基于NS模型+久期衡量利率风险（Julia）

今天刚刚交完金融风险管理的作业，我是觉得这次作业可以综合体现我上大学以来学到的东西，甚至可以代表我的水平了（让大家见笑了）。

首先，我们利用NS模型计算国债的利率期限结构。

假设顺势远期利率形式满足下式：

(1)

其中，为瞬时远期收益率，为水平因子，极限下不会衰减，代表长期成分贡献，为斜率因子，其衰减速度最快，代表短期成分贡献，为曲率因子，从0开始并衰减至0，代表中期成分贡献，为衰减率。到期收益率为瞬时远期收益率的均值，其计算公式如下： 

(2)

通过式（1）（2）计算得国债利率期限结构，作为企业债的贴现率

using XLSX#读入数据
t=XLSX.readtable("D://julia//NS.xlsx","Sheet1","A")
r=XLSX.readtable("D://julia//NS.xlsx","Sheet1","B")
t=t[1][1]
r=r[1][1]

function NS(r,t,γ)#基于NS模型的OLS回归
    x1=zeros(length(t))#用于储存x1
    for i in 1:length(t)#计算各期x1的值
        x1[i]=(γ/t[i])*(1-exp(-t[i]/γ))
    end
    x2=zeros(length(t))#用于储存x2
    for i in 1:length(t)#计算各期x2的值
        x2[i]=(γ/t[i])*(1-exp(-t[i]/γ))-exp(-t[i]/γ)
    end

    #方法一：用GLM包进行OLS回归
    #using GLM
    #using DataFrames
    #data=DataFrame(:x1=>Float64.(x1),:x2=>Float64.(x2),:r=>Float64.(r))
    #model=lm(@formula(r~x1+x2),data)

    X=[ones(length(t)) x1 x2]#将所有解释变量加常数项放入X中
    Y=r#被解释变量
    beta=(inv(X'*X))*(X')*Y

    R=zeros(length(t))#预测
    for i in 1:length(t)#依据回归后的表达式计算各剩余期限的到期收益率
        R[i]=beta[1]+beta[2]*x1[i]+beta[3]*x2[i]
    end
    error=0#计算误差
    for i in 1:length(t)
        error=error+(R[i]-r[i])^2
    end
    return error,beta,R
end

γ=0.5#设置γ的初始值
n=1#计次
Γ=zeros(60)#用于储存每次的γ
Error=zeros(60)#用于储存每次的误差，即目标函数
while γ<=30#设置γ的终值
    error,beta,R=NS(r,t,γ)
    Γ[n]=γ#储存γ
    Error[n]=error#储存误差
    γ=0.5+0.5*n#增加γ的值（这样表示可以减少计算精度问题，但不能完全避免）
    n=n+1#累计次数
end

#寻找误差最小时所对应的索引值，并输出其对应的γ
a=Error[1]
for j in 1:length(Error)
    b=Error[j]
    if b<=a
        a=b
        global m=j
    end
end
error,beta,R=NS(r,t,Γ[m])
γ=Γ[m]

接下来，用NS模型计算出来的国债利率期限结构作为贴现率，计算久期。主要是麦考利久期、FW久期、修正久期和凸度。由于学艺不精，一开始没想起来麦考利久期假设各期利率相等，写了个贴现函数，后来学FW久期刚好用上了。不过，说实话我感觉我的代码还是一如既往的效率不高。

#初步准备：
function Tiexian(r,P,t) #贴现函数，r为贴现利率序列(长度为10，从第1年到第10年)，P为现金流序列（长度为t,从对应（10-t）年份到第10年），t为时间长度
    m=zeros(t)#生成长度为t的零矩阵
    for j in 1:t#取出P中的元素，进行操作后给m赋值
        for i in (length(r)-t+1):(length(r)-t+j)#取出对应年份的利率进行操作
            P[j]=P[j]/((1+r[i])^1)#循环除以每年的利率进行贴现
        end
        m[j]=P[j]#给m赋值
    end
    return m
end
function PV(F,i,r,t)
    C=F*i#计算每月利息
    n=zeros(t)#生成长度为t的零矩阵
    for j in 1:(t-1)#计算每个月的利息
        n[j]=C
    end
    n[t]=C+F#计算最后一个月的利息
    q=Tiexian(r,n,t)#贴现得到当期债券价值
    return sum(q)
end

#麦考利久期、FW久期
function D(F,i,r,t)#麦考利久期函数，F为债券面值，i为票面利率，r为贴现利率（r相同则计算麦考利久期，r为各期的贴现率则为FW久期）,t为时间长度
    C=F*i#计算每月利息（也可以把C换成序列，这样就可以计算每月利息不同的情况，还省掉了下面的一些步骤）
    m=zeros(t)#生成长度为t的零矩阵
    for j in 1:(t-1)#计算每个月的久期分子并付给m
        m[j]=j*C#t*C,时间*利息
    end
    m[t]=t*(C+F)#最后一个单独赋值（债券本金偿还）
    p=Tiexian(r,m,t)#贴现得到麦考利久期的分子
    n=zeros(t)#生成长度为t的零矩阵
    for j in 1:(t-1)#计算每个月的利息
        n[j]=C
    end
    n[t]=C+F#计算最后一个月的利息
    q=Tiexian(r,n,t)#贴现得到当期债券价值
    return sum(p)/sum(q)#麦考利久期
end

#修正久期（麦考利久期隐含假设收益率不变，即y为定值）
function reviseD(P,i,r,t)#修正久期函数
    y=sum(r)/length(r)#计算债券到期时间内利率的均值
    reviseD=D(P,i,r,t)/(1+y)#依据修正久期和久期之间的关系
    return reviseD
end

#依据修正久期计算资产净现值变化量：
rd=reviseD(P,i,r,t)
dy=0.01
dp=-P*rd*dy

#凸度：
function C(F,i,r,t)#凸度函数
    C=F*i#计算每月利息（也可以把C换成序列，这样就可以计算每月利息不同的情况，还省掉了下面的一些步骤）
    m=zeros(t)#生成长度为t的零矩阵
    for j in 1:(t-1)#计算每个月的久期分子并付给m
        m[j]=j*(j+1)*C#t*(t+1)*C,时间*(时间+1)*利息
    end
    m[t]=t*(t+1)*(C+F)#最后一个单独赋值（债券本金偿还）
    p=Tiexian(r,m,t)#贴现得到凸度的分子
    n=zeros(t)#生成长度为t的零矩阵
    for j in 1:(t-1)#计算每个月的利息
        n[j]=C
    end
    n[t]=C+F#计算最后一个月的利息
    q=Tiexian(r,n,t)#贴现得到当期债券价值
    y=sum(r)/length(r)
    return sum(p)/(sum(q)*(1+y)^2) #凸度
end

#依据凸度计算资产净现值变化量：
c=C(P,i,r,t)
dy=0.01
dP=-P*rd*dy+0.5*P*c*(dy^2)