你要的数据,程序不是都有帮你带出?

Ex:Estimate the CIR model

EVALUATING S at FINAL PARAMETER ESTIMATES

  ------------------  GMM PARAMETER ESTIMATES  -----------------
  Parameter        Coeff    Std Err   Null   t-stat    p-val
  alpha         0.028733   0.013869   0.00     2.07   0.0383
  beta         -0.420795   0.242914   0.00    -1.73   0.0832
  sigma^2       0.006557   0.001677   0.00     3.91   0.0001

  -------------------  GMM MOMENT CONDITIONS  ------------------
                    Moment    Std Err   t-stat    p-val
     Moment 1     0.000061   0.000033     1.83   0.0679
     Moment 2    -0.000007   0.000004    -1.83   0.0680
     Moment 3     0.000029   0.000016     1.83   0.0680
     Moment 4     0.000004   0.000002     1.83   0.0680

      J-stat = 4.0487    Prob[Chi-sq.(1) > J] = 0.0442
%%%%%%%%%%%

当然你也可以自己算.

PARAMETERS Std Err
sqrt(diag(gout.bcov))
ans =

    0.0139
    0.2429
    0.0017

%%%%%%%%%%
Ex: Estimate Nine Nested Models
                   alpha    beta     Sigma^2    gamma
Full              0.0419   -0.6077    1.7788    1.5081
                ( 0.0159) ( 0.2698) ( 2.9056) ( 0.3135)

Merton            0.0051    0.0000    0.0004    0.0000
                           ( 0.0034)           ( 0.0001)         

Vasicek           0.0301   -0.4485    0.0004    0.0000
                ( 0.0141) ( 0.2468) ( 0.0001)         

CIR SR            0.0287   -0.4208    0.0066    0.5000
                ( 0.0139) ( 0.2429) ( 0.0017)         

Dothan            0.0000    0.0000    0.1061    1.0000
                                    ( 0.0219)         

GBM               0.0000    0.0772    0.1034    1.0000
                          ( 0.0603) ( 0.0224)         

Brennan-Schwarz   0.0316   -0.4527    0.1139    1.0000
                ( 0.0140) ( 0.2434) ( 0.0222)         

CIR VR            0.0000    0.0000    1.4687    1.5000
                                     ( 0.2569)         

CEV               0.0000    0.0805    0.6898    1.3579
                          ( 0.0615) ( 1.4675) ( 0.4076)

                      J_T   p-value   df    R^2_1    R^2_2
Full               0.0000              0   0.0266   0.1576
Merton             4.8451    0.0887    2   0.0000   0.0000
Vasicek            4.3807    0.0363    1   0.0145   0.0000
CIR SR             4.0487    0.0442    1   0.0128   0.0038
Dothan             7.0385    0.0707    3   0.0000   0.0282
GBM                5.7120    0.0575    2   0.0004   0.0268
Brennan-Schwarz    1.9755    0.1599    1   0.0148   0.0325
CIR VR             7.1974    0.0659    3   0.0000   0.1145
CEV                5.5478    0.0185    1   0.0005   0.0761

epoh老师，您好！
     您的估计结果用的程序和结果能粘上来，我看看吗？
这些结果不知道具体的程序和数据，所以很难看懂。
非常感谢！

哈哈zhangtao兄也有兴趣.

GMM Libraries for Matlab:

主网站目前已不提供下载

http://www.feweb.vu.nl/econometriclinks/mcliffprogs.html

但可由底下帖子免费下载

https://bbs.pinggu.org/forum.php?mod=viewthread&tid=114904&highlight=gmm

执行文件:ckls_d.m

但由于matlab版本更新

ckls_d.m请修改底下两行

% --- Estimate Nine Nested Models -------------

gmmopt.prt = 0;      -->gmmopt.prt = 1;

gmmopt.infoz.prt = 0;-->gmmopt.infoz.prt = 1;

否则fprintf会产生错误信息.

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

相同数据ckls.dat,也有S-plus code.

Modelling Financial Time Series with S-PLUS
chap21. Generalized Method of Moments
21.7.5 Interest Rate Diffusion Model
#ckls
ckls.moments <- function(parm, data = NULL,dt = 1/12) {
# parm = (alpha,beta,sigma,gamma)’
# data = [r(t+dt)-r(t),r(t)]
# dt = discretization step
e.hat = as.vector(data[,1] -(parm[1] + parm[2]*data[,2])*dt)
m2 = e.hat*as.vector(data[,2])
m3 = e.hat^2 - dt*parm[3]*parm[3]*(as.vector(data[,2])^(2*parm[4]))
m4 = m3*data[,2]
cbind(e.hat,m2,m3,m4)
}
data.ckls.ts = seriesMerge(diff(ckls.ts), tslag(ckls.ts))
colIds(data.ckls.ts)[1] = "RF.diff"
data.ckls = as.matrix(seriesData(data.ckls.ts))
summaryStats(data.ckls)
start.vals = c(0.06,-0.5,1,1)
names(start.vals) = c("alpha","beta","sigma","gamma")
gmm.ckls = GMM(start.vals,ckls.moments,ts = T,data = data.ckls,dt = 1/12)
summary(gmm.ckls)

####
#GMM estimates for the restricted models cir
cir.moments <- function(parm, data = NULL, dt = 1/12) {
# parm = (alpha,beta,sigma)’
# data = [r(t+dt)-r(t),r(t)]
# dt = discretization step
e.hat = as.vector(data[,1] -(parm[1] + parm[2]*data[,2])*dt)
m2 = e.hat*as.vector(data[,2])
m3 = e.hat^2 - dt*parm[3]*parm[3]*(as.vector(data[,2]))
m4 = m3*data[,2]
cbind(e.hat,m2,m3,m4)
}
start.vals = c(0.06,-0.5,1)
names(start.vals) = c("alpha","beta","sigma")
gmm.cir = GMM(start.vals, cir.moments, ts = T,data = data.ckls, dt = 1/12)
summary(gmm.cir,print.moments = F)

谢谢epoh老师， 不知道是否对Nowman Gaussian estimation 做CKLS 也有研究？我用matlab编写的code运行有问题，方便的话，可否指教。

Besides, 想请问你给的表格中t-stat和p-val 是通过library的哪个code得出？针对MLE method, 是否通用呢？

希望能加QQ to discuss more. Thxs!!

感谢大家的留言，以下是我写的用Nowman Gaussian Estimation做CKLS的matlab code, 但是结果运行不出来，错误显示，“??? Undefined function or method 'NEW' for input arguments of type 'char'.” 请教指点。我用的是fmincon做minimizing log likilihood value.

%Main program
function main()

opts=optimset( 'Display',             'notify', ...
                    'MaxIter',              1000,   ...
                    'MaxFunEvals',          4000,   ...
                    'LargeScale',           'off',  ...
                    'HessUpdate',           'bfgs', ...
                    'TolX',                 1e-8,   ...
                    'TolFun',               1e-11   ...
                       );   
% option settings for optimization algorithm

    load CKLSdata;
    n=length(X);
    dt=1/260;

ini_P=rand(4,1);%initial parameters [alpha,beta,sigma,gamma] as a vector.
low_P=-100*ones(4,1);%lower bounds of parameters
up_P=100*ones(4,1);%upper bounds of parameters
A2=[-1,0,0,0;0,0,0,0;0,0,-1,0;0,0,0,0];% parameter restriction for Merton
A3=[-1,0,0,0;0,1,0,0;0,0,-1,0;0,0,0,0];%parameter restriction for Vasicek
A4=[-1,0,0,0;0,1,0,0;0,0,-1,0;0,0,0,0];%parameter restriction for CIR
A5=[0,0,0,0;0,-1,0,0;0,0,-1,0;0,0,0,0];%parameter restriction for Dothan
A6=[-1,0,0,0;0,1,0,0;0,0,-1,0;0,0,0,0];%parameter restriction for Brennan and Schwartz
B=[0,0,0,0];

[P1,lnL1]=fmincon('CKLS_Calibration',ini_P,[],[],[],[],low_P,up_P,[],opts);
[P2,lnL2]=fmincon('Merton_Calibration',ini_P,A2,B,[],[],low_P,up_P,[],opts);
[P3,lnL3]=fmincon('Vasicek_Calibration',ini_P,A3,B,[],[],low_P,up_P,[],opts);
[P4,lnL4]=fmincon('CIR_Calibration',ini_P,A4,B,[],[],low_P,up_P,[],opts);
[P5,lnL5]=fmincon('Dothan_Calibration',ini_P,A5,B,[],[],low_P,up_P,[],opts);
[P6,lnL6]=fmincon('Brennan_Schwartz_Calibration',ini_P,A6,B,[],[],low_P,up_P,[],opts);

% P is parameters estimation
% lnL is minus maximized likelihood value

format long
disp(num2str([P1,lnL1],5));
disp(num2str([P2,lnL2],5));
disp(num2str([P3,lnL3],5));
disp(num2str([P4,lnL4],5));
disp(num2str([P5,lnL5],5));
disp(num2str([P6,lnL6],5));

end % end main program

function lnL1 = CKLS_Calibration(P,X) % The minus log likelihood function of CKLS model
   
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*X(i-1)^(2*P(4,1))*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1)-(P(1,1)*(exp(P(2,1)*dt)-1))/P(2,1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
    lnL1=-0.5*sum;
    lnL1=-lnL1;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end


function lnL2 = Merton_Calibration(P,X) % The minus log likelihood function of Merton model
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=(P(3,1)^2)*dt;
        ita=zeros(n,1);
        ita(i,1)=X(i)-X(i-1)-P(1,1)*dt;
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL2=-0.5*sum;
lnL2=-lnL2;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end

function lnL3 = Vasicek_Calibration(P,X)% The minus log likelihood function of Vasicek model
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1)-(P(1,1)*(exp(P(2,1)*dt)-1))/P(2,1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL3=-0.5*sum;
lnL3=-lnL3;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end



function lnL4 = CIR_Calibration(P,X) % The minus log likelihood function of CIR model
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*X(i-1)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1)-(P(1,1)*(exp(P(2,1)*dt)-1))/P(2,1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL4=-0.5*sum;
lnL4=-lnL4;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end



function lnL5 = Dothan_Calibration(P,X) % The minus log likelihood function of Dothan model
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*(X(i-1)^2)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL5=-0.5*sum;
lnL5=-lnL5;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end



function lnL6 = Brennan_Schwartz_Calibration(P,X) % The minus log likelihood function of BS model
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*(X(i-1)^2)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1)-(P(1,1)*(exp(P(2,1)*dt)-1))/P(2,1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL6=-0.5*sum;
lnL6=-lnL6;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end

先改一个供你参考,

结果跟ini_P息息相关.

%%%%%%%%%%%
load ckls.dat
data=ckls;
dt=1/12;
% option settings for optimization algorithm
opts=optimset( 'Display',             'notify', ...
                    'MaxIter',              1000,   ...
                    'MaxFunEvals',          4000,   ...
                    'LargeScale',           'off',  ...
                    'HessUpdate',           'bfgs', ...
                    'TolX',                 1e-8,   ...
                    'TolFun',               1e-11   ...
                       );   

ini_P=rand(4,1);%initial parameters [alpha,beta,sigma,gamma] as a vector.
low_P=-100*ones(4,1);%lower bounds of parameters
up_P=100*ones(4,1);%upper bounds of parameters
[P1,lnL1]=fmincon(@(P)CKLS_Calibration(P,data,dt),ini_P,[],[],[],[],low_P,up_P,[],opts)
P1 =

    0.3766
    0.0013
    0.9054
   -0.1675

lnL1 =

   -3.6264

%%%%%%%%
function lnL1 = CKLS_Calibration(P,data,dt) % The minus log likelihood function of CKLS
     n=length(data);
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3)^2)*data(i-1)^(2*P(4))*(exp(2*P(2)*dt)-1))/(2*P(2));
        ita=zeros(n,1);
        ita(i,1)=data(i)-exp(P(2)*dt)*data(i-1)-(P(1)*(exp(P(2)*dt)-1))/P(2);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
    lnL1=-0.5*sum;
    lnL1=-lnL1;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end

THXS!!!!!epoh!!!!
I LOVE YOU SOOOOOOO MUCH!!!!!!
Yeah, i got the parameters for CKLS, which turns out to be nearly the same as what i have done using Excel solver.
P1 =

   -0.1253
   -0.3716
    4.3992
    3.5434


lnL1 =

  -12.1230

I will do several others in a moment and will let you know how its going. Thxs again!!
U r such a nice person and wish you all the best!:)

epoh:

Aftering changing it like what u wrote, the first one works, but others all get stuck...
Errors show:??? Error using ==> optimfcnchk at 300
NONLCON must be a function.
Error in ==> fmincon at 382
   confcn = optimfcnchk(NONLCON,'fmincon',length(varargin),funValCheck,flags.gradconst,false,true);
Error in ==> CKLS at 25
[P1,lnL1]=fmincon(@(P)CKLS_Calibration(P,data,dt),ini_P,[],[],[],[],low_P,up_P,opts);

Hope to hear from you asap~~

Following given is my modified code:
function main()
load CKLSdata.dat
data=CKLSdata;
dt=1/260;
% option settings for optimization algorithm
opts=optimset( 'Display',             'notify', ...
                    'MaxIter',              1000,   ...
                    'MaxFunEvals',          10000,   ...
                    'LargeScale',           'off',  ...
                    'HessUpdate',           'bfgs', ...
                    'TolX',                 1e-8,   ...
                    'TolFun',               1e-11   ...
                       );   

ini_P=rand(4,1);%initial parameters [alpha,beta,sigma,gamma] as a vector.
low_P=-100*ones(4,1);%lower bounds of parameters
up_P=100*ones(4,1);%upper bounds of parameters
A2=[-1,0,0,0;0,0,0,0;0,0,-1,0;0,0,0,0];% parameter restriction for Merton
A3=[-1,0,0,0;0,1,0,0;0,0,-1,0;0,0,0,0];%parameter restriction for Vasicek
A4=[-1,0,0,0;0,1,0,0;0,0,-1,0;0,0,0,0];%parameter restriction for CIR
A5=[0,0,0,0;0,-1,0,0;0,0,-1,0;0,0,0,0];%parameter restriction for Dothan
A6=[-1,0,0,0;0,1,0,0;0,0,-1,0;0,0,0,0];%parameter restriction for Brennan and Schwartz
B=[0,0,0,0];
[P1,lnL1]=fmincon(@(P)CKLS_Calibration(P,data,dt),ini_P,[],[],[],[],low_P,up_P,opts);
[P2,lnL2]=fmincon(@(P)Merton_Calibration(P,data,dt),ini_P,[],[],[],[],low_P,up_P,[],opts);
[P3,lnL3]=fmincon(@(P)Vasicek_Calibration(P,data,dt),ini_P,[],[],[],[],low_P,up_P,[],opts);
[P4,lnL4]=fmincon(@(P)CIR_Calibration(P,data,dt),ini_P,[],[],[],[],low_P,up_P,[],opts);
[P5,lnL5]=fmincon(@(P)Dothan_Calibration(P,data,dt),ini_P,[],[],[],[],low_P,up_P,[],opts);
[P6,lnL6]=fmincon(@(P)Brennan_Schwartz_Calibration(P,data,dt),ini_P,[],[],[],[],low_P,up_P,[],opts);

end
%%%%%%%%
function lnL1 = CKLS_Calibration(P,data,dt) % The minus log likelihood function of CKLS
     n=length(data);
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3)^2)*data(i-1)^(2*P(4))*(exp(2*P(2)*dt)-1))/(2*P(2));
        ita=zeros(n,1);
        ita(i,1)=data(i)-exp(P(2)*dt)*data(i-1)-(P(1)*(exp(P(2)*dt)-1))/P(2);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
    lnL1=-0.5*sum;
    lnL1=-lnL1;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end

function lnL2 = Merton_Calibration(P,data,dt) % The minus log likelihood function of Merton model
    n=length(data);
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=(P(3,1)^2)*dt;
        ita=zeros(n,1);
        ita(i,1)=X(i)-X(i-1)-P(1,1)*dt;
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL2=-0.5*sum;
lnL2=-lnL2;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end

function lnL3 = Vasicek_Calibration(P,data,dt)% The minus log likelihood function of Vasicek model
    n=length(data);
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1)-(P(1,1)*(exp(P(2,1)*dt)-1))/P(2,1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL3=-0.5*sum;
lnL3=-lnL3;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end



function lnL4 = CIR_Calibration(P,data,dt) % The minus log likelihood function of CIR model
    n=length(data);
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*X(i-1)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1)-(P(1,1)*(exp(P(2,1)*dt)-1))/P(2,1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL4=-0.5*sum;
lnL4=-lnL4;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end



function lnL5 = Dothan_Calibration(P,data,dt) % The minus log likelihood function of Dothan model
   n=length(data);
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*(X(i-1)^2)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL5=-0.5*sum;
lnL5=-lnL5;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end



function lnL6 = Brennan_Schwartz_Calibration(P,data,dt) % The minus log likelihood function of BS model
    n=length(data);
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*(X(i-1)^2)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1)-(P(1,1)*(exp(P(2,1)*dt)-1))/P(2,1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL6=-0.5*sum;
lnL6=-lnL6;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end

更改如下:

1.我把数据名称data改回你原来的X,以免部分没更正到

   

2.似乎不需linear inequalities A*x =< B,

  我改依Nowman Ppaper page 4/13,Table 1

  照不同模型分别设置 lower &upper bound

%%%%%%%%%%%%
load ckls.dat
data=ckls(1:(end-1),3);
X=data/100;
dt=1/12;
% option settings for optimization algorithm
opts=optimset( 'Display',             'notify', ...
                    'MaxIter',              1000,   ...
                    'MaxFunEvals',          4000,   ...
                    'LargeScale',           'off',  ...
                    'HessUpdate',           'bfgs', ...
                    'TolX',                 1e-8,   ...
                    'TolFun',               1e-11   ...
                       );   

ini_P=rand(4,1);%initial parameters [alpha,beta,sigma,gamma] as a vector.
%ini_P=[0.1656 ;   0.6020;    0.2630 ;   0.65410];

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;0;-400;1];        %parameter restriction for Dothan
up_P5=[0;0;400;1];
low_P6=[-400;-400;-400;1];  %parameter restriction for Brennan and Schwartz
up_P6=[400;400;400;1];

[P1,lnL1]=fmincon(@(P)CKLS_Calibration(P,X,dt),ini_P,[],[],[],[],low_P,up_P,[],opts)
[P2,lnL2]=fmincon(@(P)Merton_Calibration(P,X,dt),ini_P,[],[],[],[],low_P2,up_P2,[],opts)
[P3,lnL3]=fmincon(@(P)Vasicek_Calibration(P,X,dt),ini_P,[],[],[],[],low_P3,up_P3,[],opts)
[P4,lnL4]=fmincon(@(P)CIR_Calibration(P,X,dt),ini_P,[],[],[],[],low_P4,up_P4,[],opts)
[P5,lnL5]=fmincon(@(P)Dothan_Calibration(P,X,dt),ini_P,[],[],[],[],low_P5,up_P5,[],opts)
[P6,lnL6]=fmincon(@(P)Brennan_Schwartz_Calibration(P,X,dt),ini_P,[],[],[],[],low_P6,up_P6,[],opts)

P1 =
    0.0346
   -0.5488
    0.6964
    3.5288
lnL1 =
  -10.6499

P2 =
   -0.0078
         0
    0.0001
         0
lnL2 =
  -10.3478

P3 =
   -0.0515
    0.5664
   -0.0001
         0
lnL3 =
  -10.3335

P4 =

   -0.0329
    0.3256
    0.0001
    0.5000
lnL4 =
  -11.9802

P5 =
         0
         0
    0.4314
    1.0000
lnL5 =

   NaN

P6 =
    0.1451
   -1.9804
   -0.0071
    1.0000
lnL6 =
   -8.8241

我把revised后的GMM library传上来.

  

内有三个文件夹 gmm,minz,Utilities

你运行时rows(ckls)出错

不是缺Utilities,就是没设好path

1.安装时,很重要!!记得点选

  file\set path\add with subfolders

2.我填加了两个文件ckls_info.m,ckls_prep.m

  方便你单独做模型时不会出错

3.执行文件:ckls_d

4.亦可单独执行如下:

5.请注意看gout=gmm(b0,gmmopt,dr,r,Z);

  gout包含了你所有需要的数据.

%  gmmout.f     function value
%  gmmout.J     chi-square stat for model fit
%  gmmout.p     p-value for model fit
%  gmmout.b     coefficient estimates
%  gmmout.se    standard errors for each parameter
%  gmmout.bcov  cov matrix of parameter estimates
%  gmmout.t     t-stats for parms = null
%  gmmout.pb    p-values for coefficients
%  gmmout.m     moments
%  gmmout.mse   standard errors of moments
%  gmmout.varm  covariance matrix of moments
%  gmmout.mt    t-stats for moments = 0
%  gmmout.mp    p-vals for moments
%  gmmout.nobs  number of observations
%  gmmout.north number of orthogonality conditions
%  gmmout.neq   number of equations
%  gmmout.nz    number of instruments
%  gmmout.nvar  number of parameters
%  gmmout.df    degrees of freedom for model

........

........

% --- Initializations ----------------------------------------------------
load ckls.dat    % Load Data
T = rows(ckls)-1;   % nb. obs.
r = ckls(1:T,3)/100;   % Rt
dr = ckls(2:end,3)/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

% --- Print Summary Stats -----------------------------------------------
in.rnames = strvcat('Mean','Std',[repmat('Rho_',6,1) int2str([1:6]')]);
in.cnames = strvcat('r(t)','r(t)-r(t-1)');
mprint1([mean([r dr]); std([r dr]); autocorr([r dr],6)],in);

% --- GMM Options -------------------------------------------------------
gmmopt.W0='Z';    % Initial Weighting Matrix
gmmopt.infoz.hess = 'gn';  % Search Direction
gmmopt.infoz.momt='ckls_dm';  % Moment Conditions
gmmopt.infoz.jake='ckls_dj';  % Jacobian of Moments
gmmopt.lags=7;       % Lags in weighting matrix
gmmopt.vname = vname;
gmmopt.null = null;
gmmopt.infoz.parms = [NaN NaN NaN NaN]'; % Special field used by ckls_dm
gmmopt.infoz
disp('These are the gmmopt and infoz structures.');
disp('Strike a key to Start Estimation.');

% --- Estimate the Full model with gmm() -------------------------------
gout=gmm(b0,gmmopt,dr,r,Z);
condvar0 = ckls_info(gout,repmat(NaN,4,1),vname,dr,r);
pause;

% --- Estimate the CIR model with gmm() --------------------------------
parms = [NaN NaN NaN .5]';  % Specify Restriction
[bin, opt] = ckls_prep(parms,b0,gmmopt,null,vname);
gout=gmm(bin,opt,dr,r,Z);
condvar = ckls_info(gout,parms,vname,dr,r);

% --- Estimate Nine Nested Models -------------------------------------
gmmopt.prt = 1;
gmmopt.infoz.prt = 1;
gmmopt.plot = 0;
Plist = strvcat('Full','Merton','Vasicek','CIR SR','Dothan','GBM',...
  'Brennan-Schwarz','CIR VR','CEV');
Parms = repmat(NaN,9,4);
out1 = Parms;
out2 = out1;
out3 = repmat(NaN,9,5);
Parms(2,[2 4]) = [0 0];  % Merton
Parms(3,4) = 0;   % Vasicek
Parms(4,4) = 0.5;  % CIR SR
Parms(5,[1 2 4]) = [0 0 1]; % Dothan
Parms(6,[1 4]) = [0 1];  % GBM
Parms(7,4) = 1;   % Brennan-Schwartz
Parms(8,[1 2 4]) = [0 0 1.5]; % CIR VR
Parms(9,1) = 0;   % CEV

for i = 1:rows(Parms)
  fprintf(1,'%s\n',Plist(i,:));
  [bin, opt] = ckls_prep(Parms(i,:)',b0,gmmopt,null,vname);
  gout=gmm(bin,opt,dr,r,Z);
  [condvar,temp] = ckls_info(gout,Parms(i,:)',vname,dr,r);
  out1(i,:) = temp(1,1:4);
  out2(i,:) = temp(2,1:4);
  out3(i,:) = [temp(1,5) temp(2,5) temp(1,6:end)];
end

rnames = Plist;
fprintf('\t\t  alpha\t\t  beta\t\t Sigma^2\t gamma\n');
for i = 1:rows(Plist)
  fprintf('%s\t',rnames(i,:));
  fprintf(' %7.4f \t %7.4f \t %7.4f \t %7.4f \n',out1(i,:));
  s = sprintf('(%7.4f)\t(%7.4f)\t(%7.4f)\t(%7.4f)\n',out2(i,:));
  s = strrep(s,'(    NaN)','         ');
  fprintf('\t\t%s\n',s);
end
in.rnames = rnames;
in.cnames = strvcat('J_T','p-value','df','R^2_1','R^2_2');
in.fmt = strvcat('%7.4f','%7.4f','%1d','%7.4f','%7.4f');
mprint1(out3,in);

epoh老师，您好！
     为什么在我的matlab7上运行会提示以下错误：
>> ckls_d
??? Error using ==> autocorr
Input 'Series' must be a vector.

Error in ==> ckls_d at 42
mprint1([mean([r dr]); std([r dr]); autocorr([r dr],6)],in);

>>

zhangtao兄:

??? Error using ==> autocorr

表示你的程序抓到garch toolbox的autocorr()

而非GMM\ Utilities\ autocorr()

请依16楼说明安装.

希望我的直觉是错的
总觉得LLF太小.
就CIR model而言
为何没见到P(4,1)=0.5
function lnL4 = CIR_Calibration(P,X,dt) % The minus log likelihood function of CIR model
    n=length(X);
    for i =2:n
        sum=0;
        v=zeros(n,1);
        v(i,1)=((P(3,1)^2)*X(i-1)*(exp(2*P(2,1)*dt)-1))/(2*P(2,1));
        ita=zeros(n,1);
        ita(i,1)=X(i)-exp(P(2,1)*dt)*X(i-1)-(P(1,1)*(exp(P(2,1)*dt)-1))/P(2,1);
        sum=sum+log(v(i,1)+(ita(i,1)^2)/v(i,1));  
    end;
lnL4=-0.5*sum;
lnL4=-lnL4;%% minimizing minus log likelihood, which is equal to maximizing log
% likelihood
end

epoh:

我的data


我觉得parameters restriction 还是有问题， 因为Vasicek, CIR, BS中应该有alpha, sigma >0, beta<0
Dothan model 原先我也是用的beta=0, 但是Prof Brigo认为应该是dr=beta*dt+sigma*dWt. 所以我是用的他建议的那个model, which is equal to GBM.
对于Dothan和Merton, 参数都应该>0.

当我把main code中的parameters restriction改成:
low_P=-100*ones(4,1);       %lower bounds of parameters
up_P=100*ones(4,1);         %upper bounds of parameters
low_P2=[0;0;0;0];     % parameter restriction for Merton
up_P2=[400;0;400;0];
low_P3=[0;-400;0;0];  %parameter restriction for Vasicek
up_P3=[400;0;400;0];
low_P4=[0;-400;0;0.5];%parameter restriction for CIR
up_P4=[400;0;400;0.5];
low_P5=[0;0;0;1];        %parameter restriction for Dothan
up_P5=[0;400;400;1];
low_P6=[0;-400;0;1];  %parameter restriction for Brennan and Schwartz
up_P6=[400;0;400;1];

带入附近中的data,结果是非常的奇怪
P1 =

   -0.1199
    0.6458
    7.6813
    3.5770


lnL1 =

  -11.7233

P2 =

         0
         0
    0.3005
         0


lnL2 =

   -3.6362

P3 =

    0.9575
         0
    0.1576
         0


lnL3 =

   NaN
P4 =

    0.9575
         0
    0.1576
    0.5000


lnL4 =

   NaN
P5 =

         0
         0
    6.6331
    1.0000


lnL5 =

   NaN
P6 =

    0.9575
         0
    0.1576
    1.0000


lnL6 =

   NaN
这个结果显然是有错误的。 你觉得问题在哪？ 难道Vasicek,CIR 那些模型也不需要给参数大于或者小于0的限制么

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:

对哦！我再把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

   

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:))

除了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, 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在你的短信息

  

%%%%%%%%%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:)

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

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不设限的话,有五个可执行

   

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];