[MATLAB] CIR利率期限结构模型的MLE估计 [推广有奖]

匿名网友

楼主

匿名网友 发表于 2015-7-8 20:37:35 |坛友微信交流群|倒序 |AI写论文

相似文件

换一批

是否 +2 论坛币

k人参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

本贴主要参考Kladıvko, K. (2007). Maximum likelihood estimation of the Cox-Ingersoll-Ross process: the Matlab implementation. Technical Computing Prague.

CIR 利率期限结构模型

\[d{{r}_{t}}=\alpha (\mu -{{r}_{t}})dt+\sqrt{{{r}_{t}}}\sigma d{{W}_{t}}\]

论文及数据集Pribor3M见附件：

clc
clear all
close all
load Pribor3M
Model.Data = Pribor3M; %此处即是文中Model.Data
Model.TimeStep = 1/250; % recommended: 1/250 for daily data, 1/12 for monthly data, etc
Model.Disp = 'y'; % 'y' | 'n' (default: y)
Model.MatlabDisp = 'iter'; % 'off'|'iter'|'notify'|'final' (default: off)
Model.Method = 'besseli'; % 'besseli' | 'ncx2pdf' (default: besseli)
Results = CIRestimation(Model);
% 使用OLS估计alpha mu sigma
Nobs = length(Model.Data);
x = Model.Data(1:end-1);
dx = diff(Model.Data);
dx = dx./x.^0.5;
regressors = [Model.TimeStep./x.^0.5, Model.TimeStep*x.^0.5];
drift = regressors\dx;
res = regressors*drift - dx;
alpha = -drift(2);
mu = -drift(1)/drift(2);
sigma = sqrt(var(res, 1)/Model.TimeStep);
InitialParams = [alpha mu sigma];
if ~isfield(Model, 'Disp'), Model.Disp = 'y'; end;
if strcmp(Model.Disp, 'y')
fprintf('\n initial alpha = %+3.6f\n initial mu = %+3.6f\n initial sigma = %+3.6f\n', alpha, mu, sigma);
end
% Optimization using fminsearch 优化过程
if ~isfield(Model, 'MatlabDisp'), Model.MatlabDisp = 'off'; end;
options = optimset('LargeScale', 'off', 'MaxIter', 300, 'MaxFunEvals', 300, 'Display', Model.MatlabDisp, 'TolFun', 1e-4, 'TolX', 1e-4, 'TolCon', 1e-4);
if ~isfield(Model, 'Method'), Model.Method = 'besseli'; end;
if strcmp(Model.Method, 'ncx2pdf')
[Params, Fval, Exitflag] = fminsearch(@(Params) CIRobjective2(Params, Model), InitialParams, options);
else
[Params, Fval, Exitflag] = fminsearch(@(Params) CIRobjective1(Params, Model), InitialParams, options);
end
Results.Params = Params;
Results.Fval = -Fval/Nobs;
Results.Exitflag = Exitflag;
if strcmp(Model.Disp, 'y')
fprintf('\n alpha = %+3.6f\n mu = %+3.6f\n sigma = %+3.6f\n', Params(1), Params(2), Params(3));
fprintf(' log-likelihood = %+3.6f\n', -Fval/Nobs);
end
%最大似然估计方法之一：贝塞尔函数
Data = Model.Data;
DataF = Data(2:end);
DataL = Data(1:end-1);
Nobs = length(Data);
TimeStep = Model.TimeStep;
alpha = Params(1);
mu = Params(2);
sigma = Params(3);
c = 2*alpha/(sigma^2*(1-exp(-alpha*TimeStep)));
q = 2*alpha*mu/sigma^2-1;
u = c*exp(-alpha*TimeStep)*DataL;
v = c*DataF;
z = 2*sqrt(u.*v);
bf = besseli(q,z,1);
lnL= -(Nobs-1)*log(c) + sum(u + v - 0.5*q*log(v./u) - log(bf) - z);
%最大似然估计之二：非中心卡方分布函数
Data = Model.Data;
DataF = Data(2:end);
DataL = Data(1:end-1);
TimeStep = Model.TimeStep;
alpha = Params(1);
mu = Params(2);
sigma = Params(3);
c = 2*alpha/(sigma^2*(1-exp(-alpha*TimeStep)));
q = 2*alpha*mu/sigma^2-1;
u = c*exp(-alpha*TimeStep)*DataL;
v = c*DataF;
nc = 2*u;
df = 2*q+2;
s = 2*v;
gpdf = ncx2pdf(s, df, nc);
ppdf = 2*c*gpdf;
lnL = sum(-log(ppdf));