哪位大神能帮着看一下对不对
用最大似然法求CIR参数
Model.dett = 1/360;
Nobs = length (Model.Data);
r = Model.Data(1:Nobs-1);
dr = diff(Model.Data)./r.^0.5;
reg = [Model.dett./r.^0.5, Model.dett.*r.^0.5];
drift = reg\dr;
res = reg*drift-dr;
alpha = -drift(2);
mu = -drift(1)/drift(2);
sigma = sqrt(var(res, 1)/ Model.dett);
Params = [alpha mu sigma];

Data = Model.Data;
Nobs = length(Data);
DataA = Data(2:Nobs);
DataB = Data(1:Nobs-1);

dett=Model.dett;
alpha = Params(1);
mu = Params (2);
sigma = Params(3);
c=2*alpha/(sigma^2*(1-exp(-alpha*dett)));
q=2*alpha*mu/sigma^2-1;
u=c*exp(-alpha*dett)*DataB;
v=c*DataA;

nc=2*u;
df=2*q+2;
s=2*v;

gpdf = ncx2pdf(s, df, nc);
ppdf = 2*c*gpdf;
lnL = sum(-log(ppdf));

Nobs = length (Model.Data);
z=2*sqrt(u.*v);
bf=besseli(q,z,1);
lnM=-(Nobs-1)*log(c)+sum(u+v-0.5*q*log(v./u)-log(bf)-z);

InitialParams = Params;
if ~isfield(Model,'MatlabDisp'), Model.MatlabDisp = 'off';end;
options = optimset('LargeScale','off', 'MaxIter', 300, 'MaxFunEvals',300, 'Display', Model.MatlabDisp, 'TolFun', 1e-6, 'TolX', 1e-4, 'TolCon', 1e-6);
if ~isfield(Model,'Method'), Model.Method = 'besseli';end;
if strcmp(Model.Method, 'ncx2pdf')
    [Params, Fval, Exitflag] = fminsearch(@(Params) lnL, InitialParams,options);
else
    [Params, Fval, Exitflag] = fminsearch(@(Params) lnM, InitialParams,options);
end
Objvalue = [Fval, Exitflag];
outParams = [Params(1),Params(2),Params(3)];