跪求CKLS GMM法参数估计T_T

epoh 发表于 2011-9-1 20:37
希望我的直觉是错的
总觉得LLF太小.
就CIR model而言

hey, epoh!

sorry that the data i attached before was wrong. I sent you a new one to your yahoo email address.

Plz check that:))

刚看了下Interest Rate Models_ Theory and Practice

发现跟你说的限制都不同

或许是公式不同

麻烦你看一下

就以Dothan而言:

3.2.2 The Dothan Model
dr(t) = ar(t)dt + σr(t)dW(t), (3.17)
where a is a real constant
low_P5=[0;-400;0;1];        %parameter restriction for Dothan
up_P5=[0;400;400;1];
[P5,lnL5]=fmincon(@(P)Dothan_Calibration(P,X,dt),ini_P,[],[],[],[],low_P5,up_P5,[],opts)
P5 =
         0
   -3.1342
    0.0006
    1.0000

lnL5 =

  -13.2832

epoh 发表于 2011-9-1 22:16
刚看了下Interest Rate Models_ Theory and Practice发现跟你说的限制都不同或许是公式不同麻烦你看一下就以 ...

epoh:

对哦！我再把restriction好好研究一下。 多谢！ data 发给你邮箱了。如果方便的话可以麻烦你用GMM 和MLE运行一下么？

我对matlab的data 调用不大懂，当我改成自己的data, in_sample.dat(copied to work space and saved as dat file) , 把程序改成如下，但是调用不出来， 显示error??? Error using ==> load
Unable to read file in_sample.dat: No such file or directory.
Error in ==> ckls_d at 29


load in_sample.dat                        % Load Data
load in_sample.dat                        % Load Data
T = rows(in_sample)-1;                        % nb. obs.
r = in_sample(1:T)/100;                        % Rt
dr = in_sample(2:end)/100-r;                % Rt+1 - Rt
Z = [ones(T,1) r];                        % Instruments
b0=[mean(r); 0; 12*var(dr)/var(r); .5];        % Parameter Starting Values
null = [0 0 0 .5]';                        % Null hypothesis
vname = strvcat('alpha','beta','sigma^2','gamma');        % Names

in_sample.dat,out_sample.dat

我已帮你转换为matlab的数据型态

in_sample.mat,out_sample.mat

    sample.mat.rar (4.1 KB)

2011-9-1 23:12:40 上传

MLE的部分明日再给你.

%%%%%%%%%%%

load in_sample
T = rows(in_sample)-1;                        % nb. obs.
r = in_sample(1:T)/100;                       % Rt
dr = in_sample(2:end)/100-r;                  % Rt+1 - Rt
Z = [ones(T,1) r];                            % Instruments
b0=[mean(r); 0; 12*var(dr)/var(r); .5];        % Parameter Starting Values
null = [0 0 0 .5]';                            % Null hypothesis
vname = strvcat('alpha','beta','sigma^2','gamma');        % Names
......
......
result:
% --- Estimate the Full model with gmm()
  ------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter         Coeff      Std Err   Null   t-stat    p-val
  alpha          0.014720     0.007334   0.00     2.01   0.0447
  beta          -0.713333     0.373422   0.00    -1.91   0.0561
  sigma^2      729.803367   610.598439   0.00     1.20   0.2320
  gamma          2.222879     0.128451   0.50    13.41   0.0000
% --- Estimate the CIR model with gmm()
------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter        Coeff    Std Err   Null   t-stat    p-val
  alpha         0.006961   0.004417   0.00     1.58   0.1151
  beta         -0.336910   0.241542   0.00    -1.39   0.1631
  sigma^2       0.002539   0.000537   0.00     4.73   0.0000

epoh 发表于 2011-9-1 23:14
in_sample.dat,out_sample.dat我已帮你转换为matlab的数据型态in_sample.mat,out_sample.mat   MLE的部分明 ...

epoh:))

除了ckls_d.m第一段要改， 其他的呢？ 因为他原来用的数据是monthly的（主文件中的解释）
% A discrete time (monthly) approximation is
%    r(t) - r(t-1) = kappa[theta - r(t-1)]*[t - (t-1)] + e(t)
% or
%     y(t) = (a + b*r(t-1))/12 + e(t)
% so
%     e1(t) = y(t) - (a + b*r(t-1))/12
%     e2(t) = e1(t)^2 - sigma^2 r^(2*gamma)/12

我是daily set的data, 不像MLE我只需要改dt, 用GMM的话，需要改其他的么？

epoh 发表于 2011-9-1 23:14
in_sample.dat,out_sample.dat我已帮你转换为matlab的数据型态in_sample.mat,out_sample.mat   MLE的部分明 ...

epoh, i am coming to bother you again, hah, sorry.

GMM:
我将主文件和ckls.info, ckls.dj,ckls.dm 中所有的12替换成了260，删掉了主程序中plot的那段（因为运行时那块有错）
结果：
                  alpha                  beta                 Sigma^2         gamma
Full                     0.0147          -0.7133          729.8034           2.2229
                ( 0.0073)        ( 0.3734)        (610.5984)        ( 0.1285)

Merton                   0.0008           0.0000           0.0001           0.0000
                ( 0.0004)                         ( 0.0000)                

Vasicek                  0.0070          -0.3384           0.0000           0.0000
                ( 0.0045)        ( 0.2444)        ( 0.0000)                

CIR SR                   0.0070          -0.3369           0.0025           0.5000
                ( 0.0044)        ( 0.2415)        ( 0.0005)                

Dothan                   0.0000           0.0000           0.1151           1.0000
                                                  ( 0.0222)                

GBM                      0.0000           0.0372           0.1221           1.0000
                                 ( 0.0245)        ( 0.0230)                

Brennan-Schwarz          0.0062          -0.2975           0.1156           1.0000
                ( 0.0044)        ( 0.2382)        ( 0.0226)                

CIR VR                   0.0000           0.0000           4.5565           1.5000
                                                  ( 0.7882)                

CEV                      0.0000           0.0322          240.0645           2.0626
                                 ( 0.0309)        (409.4776)        ( 0.2534)

                      J_T   p-value   df    R^2_1    R^2_2
Full               0.0000              0   0.0001   0.0006
Merton             2.7400    0.2541    2   0.0000   0.0000
Vasicek            1.9409    0.1636    1   0.0000   0.0000
CIR SR             1.8957    0.1686    1   0.0000   0.0000
Dothan             5.3158    0.1501    3   0.0000   0.0000
GBM                3.1978    0.2021    2   0.0000   0.0000
Brennan-Schwarz    2.1732    0.1404    1   0.0000   0.0000
CIR VR             4.2881    0.2320    3   0.0000   0.0000
CEV                3.1268    0.0770    1   0.0000   0.0003
此结果 比 不替换12的结果 中的gamma 数值更加接近我excel中的数值， 但也有问题： Vasicek模型中sigmma=0, 其他值看上去都还比较真实。

MLE:  出来的数据还是明显不对，可能是我程序的问题。 我把所用的计算lnL的公式发你yahoo邮箱了~
你有没有其他用MLE做CKLS的code, 我这个程序问题不知道是错在哪里了哎。

期待你的回复。。。。。

先说一下,你的论文内容写的太好了!

你的数据GMM都OK,

我重新写了个6个模型的 GMM code.txt.

%%%%%%%%

MLE CIR model,我找了篇文献

MAXIMUM LIKELIHOOD ESTIMATION OF THE COX-INGERSOLL-ROSS PROCESS: THE MATLAB IMPLEMENTATION

log-likelihood function的确和你的有所不同

我用你的数据做的结果如下:

lnL4 = -8.1672e+003

   parameters      se           t-value        p-value

    10.5529        0.0049       2161.3001       0.0000

     0.0206        0.0006       37.0483         0.0000

     0.2459        0.0001       3990.1709       0.0000

GMM:

Long-run mean,    theta = 2.0662%

Speed of adj,     kappa = 0.3369

Volatility parm,  sigma = 0.0504

MLE  VASICEK model,我亦找了篇文献

VASICEK INTEREST RATE MODEL:PARAMETER ESTIMATION,

EVOLUTION OF THE SHORT-TERM INTEREST RATE AND TERM STRUCTURE

这是用EXCEL做的

你可参考公式(14),(15),(17)

有了两篇文献可参考,

MLE的部分你应该很快就可修改完成

有关的程序文献一并传上

但因MLE部分可能还要修改

未臻成熟,为免误导,我设了pw.

pw在你的短信息

   GMM_MLE.rar (602.61 KB)

2011-9-2 15:30:05 上传

%%%%%%%%%Result:
%%%%%%%%%Full
------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter         Coeff      Std Err   Null   t-stat    p-val
  alpha          0.014720     0.007334   0.00     2.01   0.0447
  beta          -0.713333     0.373422   0.00    -1.91   0.0561
  sigma^2      729.803367   610.598439   0.00     1.20   0.2320
  gamma          2.222879     0.128451   0.50    13.41   0.0000

Long-run mean,    theta = 2.0636%
Speed of adj,       kappa = 0.7133
Volatility parm,    sigma = 27.0149
Cond. Vol. parm,  gamma = 2.2229
Average Cond Volatility = 0.0377%
R^2 (yld change)   = 0.0001

%%%%%%%%Merton
------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter       Coeff    Std Err   Null   t-stat    p-val
  alpha          0.000790   0.000425   0.00     1.86   0.0633
  sigma^2      0.000051   0.000011   0.00     4.60   0.0000

Constraints: beta= 0.0000 gamma= 0.0000
Long-run mean,    theta =   -Inf%
Speed of adj,       kappa = -0.0000
Volatility parm,    sigma = 0.0072
Cond. Vol. parm,  gamma = 0.0000
Average Cond Volatility = 0.0445%
R^2 (yld change)   = 0.0000

%%%%%%%Vasicek
------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter        Coeff    Std Err   Null   t-stat    p-val
  alpha         0.006980   0.004472   0.00     1.56   0.1186
  beta         -0.338391   0.244355   0.00    -1.38   0.1661
  sigma^2       0.000049   0.000011   0.00     4.43   0.0000

Constraints: gamma= 0.0000
Long-run mean,    theta = 2.0626%
Speed of adj,       kappa = 0.3384
Volatility parm,    sigma = 0.0070
Cond. Vol. parm,  gamma = 0.0000
Average Cond Volatility = 0.0434%
R^2 (yld change)   = 0.0000

%%%%%%%CIR
------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter        Coeff    Std Err   Null   t-stat    p-val
  alpha         0.006961   0.004417   0.00     1.58   0.1151
  beta         -0.336910   0.241542   0.00    -1.39   0.1631
  sigma^2       0.002539   0.000537   0.00     4.73   0.0000

Constraints: gamma= 0.5000
Long-run mean,    theta = 2.0662%
Speed of adj,       kappa = 0.3369
Volatility parm,    sigma = 0.0504
Cond. Vol. parm,  gamma = 0.5000
Average Cond Volatility = 0.0441%
R^2 (yld change)   = 0.0000

%%%%%%%Dothan
------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter       Coeff    Std Err   Null   t-stat    p-val
  sigma^2      0.115149   0.022191   0.00     5.19   0.0000

Constraints: alpha= 0.0000 beta= 0.0000 gamma= 1.0000
Long-run mean,    theta =    NaN%
Speed of adj,       kappa = -0.0000
Volatility parm,    sigma = 0.3393
Cond. Vol. parm,  gamma = 1.0000
Average Cond Volatility = 0.0435%
R^2 (yld change)   = 0.0000

%%%%%%%Brennan and Schwartz
------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter        Coeff    Std Err   Null   t-stat    p-val
  alpha         0.006171   0.004363   0.00     1.41   0.1573
  beta         -0.297537   0.238163   0.00    -1.25   0.2116
  sigma^2       0.115576   0.022563   0.00     5.12   0.0000

Constraints: gamma= 1.0000
Long-run mean,    theta = 2.0739%
Speed of adj,       kappa = 0.2975
Volatility parm,    sigma = 0.3400
Cond. Vol. parm,  gamma = 1.0000
Average Cond Volatility = 0.0436%
R^2 (yld change)   = 0.0000

epoh 发表于 2011-9-2 15:32
先说一下,你的论文内容写的太好了!你的数据GMM都OK,我重新写了个6个模型的 GMM code.txt.GMM: alpha,beta,  ...

epoh:)

谢谢你的评价，导师意见回复来了，给了一堆要修改的地方。。。头大。。。。

GMM 运行的结果和我之前替换12成260的结果一致， 但是Vasicek model的sigma还是等于0. 这个不太正常吧？！

谢谢你给的Vasicek和CIR的两篇文献， 我刚看了一下，他使用另外一种方法得到的lnL, 不是Nowman Gauusian. 所以我觉得可能不太适合。。。另外一个问题是，如果用MLE, 不知道最后怎样也得到J_T   p-value    R^2_1    R^2_2 的数值， 我对这块statistics的不是很了解， 我在文章中该如何解释这些呢？公式之类的，因为导师的意见是 当看官什么都不懂。。。

谢！

Vasicek Model,sigma^2  不是0

------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter        Coeff    Std Err   Null   t-stat    p-val
  alpha            0.006980   0.004472   0.00     1.56   0.1186
  beta            -0.338391   0.244355   0.00    -1.38   0.1661
  sigma^2       0.000049   0.000011   0.00     4.43   0.0000

Constraints: gamma= 0.0000
Long-run mean,    theta = 2.0626%
Speed of adj,       kappa = 0.3384
Volatility parm,    sigma = 0.0070
Cond. Vol. parm,  gamma = 0.0000
Average Cond Volatility = 0.0434%
R^2 (yld change)   = 0.0000
%%%%%%%%%%%%%%%%%%

就 CIR model而言

J-stat = 1.8957    Prob[Chi-sq.(1) > J] = 0.1686

算法如下:

nobs = rows(dr); %1715

stat.f=0.00110536

J = nobs*stat.f  %1.8957

然stat.f 是由function minz()计算求得

所以MLE无法计算你要的数据.

%%%%%%%%%%%%%%%%%%%

说的也是,你是要Nowman Gauusian

那MLE的部分就再依Restriction及initial values

研究一下.

检查一下你的log likelihood function

的确有些BUG,

sum=0;要放在loop前面

否则i=1,归零

否则i=2,又归零

sum=0;

for i =2:n

   ......     

end;

不是 sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));

应该--> sum=sum+log(v(i,1))+(ita(i,1)^2)/v(i,1);

有些程序有误,我有修改.

也麻烦你顺便检查一次

六个model不设限的话,有五个可执行

    REVISED_LLF.rar (2.34 KB)

2011-9-2 22:56:11 上传

ini_P=[0.421761282626275;0.915735525189067;0.792207329559554;0.959492426392903;]
low_P=-100*ones(4,1);       %lower bounds of parameters
up_P=100*ones(4,1);         %upper bounds of parameters
low_P2=[-400;0;-400;0];     % parameter restriction for Merton
up_P2=[400;0;400;0];
low_P3=[-400;-400;-400;0];  %parameter restriction for Vasicek
up_P3=[400;400;400;0];
low_P4=[-400;-400;-400;0.5];%parameter restriction for CIR
up_P4=[400;400;400;0.5];
low_P5=[0;-400;-400;1];     %parameter restriction for Dothan
up_P5=[0;400;400;1];
low_P6=[-400;-400;-400;1];  %parameter restriction for Brennan and Schwartz
up_P6=[400;400;400;1];