[matlab]用lsqcurvefit或lsqnonlin实现NSS利率期限模型-经管之家官网！

献给梁博和各位：
Matlab Fixed Income Toolbox提供了fitSvensson和fitNelsonSiegel等过程，但是这个过程并不适合我们，主要以下几个因素：
1. 需要输入完整的债券信息（起息日，到期日，息票利率），比较麻烦。并且不同的债券有不同的付息频率，而系统提供的过程只能选择整体债券的compounding频率，不符合我国国情。
2. 该过程算zero rate比较合适，而国内即期利率的运用还不是很发达呀。故它提供的过程可能一般性不高。
其实只要使用matlab提供的原始的lsqcurvefit和lsqnonlin即可算出NSS模型，拥有每个债券的Yield(zero, par, coupon, ytm)和time to maturity就可以计算各种收益率曲线了。下面的程序计算的是“到期收益率曲线”：
用lsqcurvefit实现：
第一个文件NSScurvefit.m

1. opt = optimset('lsqcurvefit');
   
2. opt.Display = 'final';
   
3. opt.MaxFunEvals=20000;
   
4. opt.MaxIter=20000;
   
5. opt.TolFun=1e-140;
   
6. opt.TolX=1e-140;
   
7. 
   
8. % x=债券还有多久到期（年）；y=对应债券到期收益率。
   
9. % 所有数据均来自于上海证券交易所上市的国债。
   
10. x = [0.0274 0.0356 0.189 0.2 0.2411 0.3808 0.3945 0.463 0.4767 0.5096 0.5151 0.5534 0.5644 0.7014 0.7068 0.8164 0.8603 0.874 0.9753 1.0521 1.0685 1.326 1.3781 1.4712 1.5644 1.589 1.7808 1.8329 1.8767 1.8795 2.011 2.074 2.1315 2.2685 2.3562 2.4192 2.5616 2.5863 2.5918 2.7644 2.811 2.8986 3.0712 3.0904 3.2877 3.3288 3.4356 3.5479 3.5507 3.6082 3.6849 3.7973 3.8192 3.9726 4.0658 4.0877 4.2849 4.2986 4.4055 4.4603 4.6082 4.6521 4.7479 4.8055 4.9068 4.9726 4.9781 5.1315 5.3616 5.3918 5.4767 5.5726 5.6521 5.726 5.8027 5.8904 5.9562 6.0712 6.2438 6.3397 6.3863 6.6466 6.6849 6.8959 6.9534 7.126 7.3644 7.5178 7.6329 7.7288 7.8822 8.0164 8.2658 8.4027 8.537 8.6521 8.7671 8.8822 8.9973 9.1315 9.2274 9.3808 9.6301 9.7534 9.8027 10.0466 11.3342 11.4658 12.0795 12.6192 13.5452 14.6603 15.8 16.789 17.3151 17.8329 18.4658 18.8493 19.5014 25.5589 26.5342 26.9945 27.4548 27.9726 28.3507 28.6575 28.7616 28.811 29.126 29.337 29.663 48.1123 48.5945 49.0822 49.6];
    
11. y = [2.6744 2.6835 2.8781 2.2192 2.9456 3.0939 3.1034 1.9007 3.142 3.1447 3.1448 3.1451 2.9415 3.1461 3.1461 2.8037 3.1465 3.1465 3.1465 3.1477 2.6143 2.9204 3.2046 3.2315 3.2611 3.2691 3.3328 3.3495 3.363 2.9622 3.4001 3.4145 3.4257 3.4461 3.4554 3.4606 3.4693 3.3315 3.4708 3.4788 3.4813 3.4867 3.5015 3.5031 3.5164 3.1026 3.5235 3.5278 3.5279 3.5299 3.5324 3.5363 3.5371 3.5435 3.5481 3.5493 3.5603 3.5611 3.5676 3.5708 3.58 3.5828 3.5888 3.5924 3.5987 3.6027 3.6031 3.614 3.6375 3.6411 3.6504 3.6608 3.6688 3.6755 3.6814 3.3918 3.6889 3.6906 3.6919 3.6925 3.6927 3.6934 3.6934 3.6936 3.6936 3.696 3.7122 3.7269 3.7382 3.7476 3.7602 3.7678 3.7766 3.7805 3.7836 3.7857 3.7875 3.7889 3.7901 3.7911 3.7917 3.9571 3.7937 3.5953 3.7946 3.7963 3.8321 3.6299 3.8695 3.9011 3.8749 4.0265 4.0778 4.1048 4.1146 4.1226 4.1311 4.1363 4.1465 4.2247 4.2304 4.2329 4.2355 4.2386 4.241 4.2431 4.2439 4.2442 4.2466 4.2482 4.2511 4.437 4.4413 4.4457 4.4505];
    
12. 
    
13. 
    
14. % lsqcurvefit不需要给出对应的方差值，只要给出y值，matlab会自动帮你计算Least square。
    
15. % lsqnonlin返回的必须是yield本身，因为lsqnonlin求的是函数值的极值。
    
16. % 当然，使用lsqnonlin可以比lsqcurvefit多一个残差选项选择，可以针对不同的期限的残差选取不同的权重。
    
17. p=lsqcurvefit(@NSScurvefit2,[0 0 0 0 4 3],x,y,[0 -inf -inf -inf 0 0],[inf inf inf inf inf inf],opt);
    
18. 
    
19. % newx是新的横坐标。
    
20. newx=0.1:0.1:50;
    
21. newy=NSScurvefit2(p,newx);
    
22. plot(newx, newy);
    
23. hold all;
    
24. plot(x,y);

复制代码第二个文件NSScurvefit2.m

1. function [f]=NSScurvefit2(c,x)
   
2. % 由于输入的x是行向量，在这里必须首先把其转换成列向量。
   
3. % 程序中的"./"千万不能忽略，曾经因为只写了一个"/"而迷糊了很久……
   
4. %原NSS公式中beta0=c(1), beta1=c(2), beta2=c(3), beta4=c(4), tau1=(c5), tau2=c(6)
   
5. y = c(1)+c(2)*(1-exp(-x(:)/c(5)))./(x(:)/c(5))...
   
6. +c(3)*((1-exp(-x(:)/c(5)))./(x(:)/c(5))-exp(-x(:)/c(5)))...
   
7. +c(4)*((1-exp(-x(:)/c(6)))./(x(:)/c(6))-exp(-x(:)/c(6)));
   
8. f=y';%最后计算完成后，再将y换成行向量。
   
9. end

复制代码
用lsqnonlin实现（微小差别）：
第一个文件NSSnonlin.m:

1. opt = optimset('lsqcurvefit');
   
2. opt.Display = 'final';
   
3. opt.MaxFunEvals=20000;
   
4. opt.MaxIter=20000;
   
5. opt.TolFun=1e-140;
   
6. opt.TolX=1e-140;
   
7. 
   
8. % lsqcurvefit不需要给出对应的方差值，只要给出y值，matlab会自动帮你计算Least square。
   
9. % lsqnonlin返回的必须是yield本身，lsqnonlin求的是函数值的极值。
   
10. % 当然，使用lsqnonlin可以比lsqcurvefit多一个残差选项选择，可以针对不同的期限选取不同的权重。
    
11. p=lsqnonlin(@NSSnonlin2,[0 0 0 0 4 3],[0 -inf -inf -inf 0 0],[inf inf inf inf inf inf],opt);
    
12. 
    
13. % newx是新的横坐标。
    
14. newx=0.1:0.1:50;
    
15. newy=NSScurvefit2(p,newx);
    
16. plot(newx, newy);
    
17. hold all;
    
18. plot(x,y);

复制代码第二个文件NSSnonlin2.m:

1. function [f]=NSSnonlin2(c)
   
2. % x=债券还有多久到期；y=对应债券到期收益率。
   
3. % 所有数据均来自于上海证券交易所上市的国债。
   
4. x = [0.0274 0.0356 0.189 0.2 0.2411 0.3808 0.3945 0.463 0.4767 0.5096 0.5151 0.5534 0.5644 0.7014 0.7068 0.8164 0.8603 0.874 0.9753 1.0521 1.0685 1.326 1.3781 1.4712 1.5644 1.589 1.7808 1.8329 1.8767 1.8795 2.011 2.074 2.1315 2.2685 2.3562 2.4192 2.5616 2.5863 2.5918 2.7644 2.811 2.8986 3.0712 3.0904 3.2877 3.3288 3.4356 3.5479 3.5507 3.6082 3.6849 3.7973 3.8192 3.9726 4.0658 4.0877 4.2849 4.2986 4.4055 4.4603 4.6082 4.6521 4.7479 4.8055 4.9068 4.9726 4.9781 5.1315 5.3616 5.3918 5.4767 5.5726 5.6521 5.726 5.8027 5.8904 5.9562 6.0712 6.2438 6.3397 6.3863 6.6466 6.6849 6.8959 6.9534 7.126 7.3644 7.5178 7.6329 7.7288 7.8822 8.0164 8.2658 8.4027 8.537 8.6521 8.7671 8.8822 8.9973 9.1315 9.2274 9.3808 9.6301 9.7534 9.8027 10.0466 11.3342 11.4658 12.0795 12.6192 13.5452 14.6603 15.8 16.789 17.3151 17.8329 18.4658 18.8493 19.5014 25.5589 26.5342 26.9945 27.4548 27.9726 28.3507 28.6575 28.7616 28.811 29.126 29.337 29.663 48.1123 48.5945 49.0822 49.6];
   
5. y = [2.6744 2.6835 2.8781 2.2192 2.9456 3.0939 3.1034 1.9007 3.142 3.1447 3.1448 3.1451 2.9415 3.1461 3.1461 2.8037 3.1465 3.1465 3.1465 3.1477 2.6143 2.9204 3.2046 3.2315 3.2611 3.2691 3.3328 3.3495 3.363 2.9622 3.4001 3.4145 3.4257 3.4461 3.4554 3.4606 3.4693 3.3315 3.4708 3.4788 3.4813 3.4867 3.5015 3.5031 3.5164 3.1026 3.5235 3.5278 3.5279 3.5299 3.5324 3.5363 3.5371 3.5435 3.5481 3.5493 3.5603 3.5611 3.5676 3.5708 3.58 3.5828 3.5888 3.5924 3.5987 3.6027 3.6031 3.614 3.6375 3.6411 3.6504 3.6608 3.6688 3.6755 3.6814 3.3918 3.6889 3.6906 3.6919 3.6925 3.6927 3.6934 3.6934 3.6936 3.6936 3.696 3.7122 3.7269 3.7382 3.7476 3.7602 3.7678 3.7766 3.7805 3.7836 3.7857 3.7875 3.7889 3.7901 3.7911 3.7917 3.9571 3.7937 3.5953 3.7946 3.7963 3.8321 3.6299 3.8695 3.9011 3.8749 4.0265 4.0778 4.1048 4.1146 4.1226 4.1311 4.1363 4.1465 4.2247 4.2304 4.2329 4.2355 4.2386 4.241 4.2431 4.2439 4.2442 4.2466 4.2482 4.2511 4.437 4.4413 4.4457 4.4505];
   
6. 
   
7. % 由于输入的x是行向量，在这里必须首先把其转换成列向量。
   
8. % 程序中的"./"千万不能忽略，曾经因为只写了一个"/"而迷糊了很久……
   
9. %原NSS公式中beta0=c(1), beta1=c(2), beta2=c(3), beta4=c(4), tau1=(c5), tau2=c(6)
   
10. t =y'-( c(1)+c(2)*(1-exp(-x(:)/c(5)))./(x(:)/c(5))...
    
11. +c(3)*((1-exp(-x(:)/c(5)))./(x(:)/c(5))-exp(-x(:)/c(5)))...
    
12. +c(4)*((1-exp(-x(:)/c(6)))./(x(:)/c(6))-exp(-x(:)/c(6))));
    
13. 
    
14. %计算残差，可以在前面那个公式里赋权重。
    
15. f=t*t';
    
16. 
    
17. end

复制代码最后的图片都是一样的：
当然，Term structure的制作过程很复杂，包括挑选债券，剔除掉交易不频繁的；对于某些银行间债券选取双边报价价差较小的作为定价基础；对于某些期限采取其他定价方式等等。
中债登给出的收益率曲线是三次Hermite插值法，并且是分段制作的，很麻烦。有兴趣的人可以继续研究一下。