format long;

total=2:168;

totalminus1=1:167;

T=168;

%T=168，total表示2到168，totalminus1表示1到167，一会需要用到。

%注意要估计的向量为theta，theta=(c,phi,sigma2)',需要写出Likelyhood

function

%在写LF之前，需要给定c,phi,sigma的initial guess。故inital

guess的theta0如下

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

c=0.5;phi=0.5;sigma2=0.5;

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

diffc=0.01;

diffphi=0.01;

diffsigma2=0.01;

tolerent=0.0001;

if diffc

c

phi

sigma2

break

else

theta(1,1)=c;

theta(2,1)=phi;

theta(3,1)=sigma2;

%下面写出Likelyhood function,注意容易写错，注意检查。

%以下是根据c，phi，sigma2来计算ltheta的过程。

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

ltheta=-0.5*log(2*pi)-0.5*log(sigma2/(1-phi^2))-((yt(1)-c/(1-phi))^2)/(2*sigma2/(1-phi^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigma2)-(1/(2*sigma2))*((yt(total)'*yt(total))-2*c*(ones(167,1)'*yt(total))-2*phi*(yt(totalminus1)'*yt(total))+((T-1)*(c^2))+2*phi*c*(ones(167,1)'*yt(totalminus1))+(phi^2)*(yt(totalminus1)'*yt(totalminus1)));

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

%%%%%%%%%%计算ltheta对c的导数%%%%%%%%%%%%

cplusdiff=c+0.000001;

lthetacplussdiff=-0.5*log(2*pi)-0.5*log(sigma2/(1-phi^2))-((yt(1)-cplusdiff/(1-phi))^2)/(2*sigma2/(1-phi^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigma2)-(1/(2*sigma2))*((yt(total)'*yt(total))-2*cplusdiff*(ones(167,1)'*yt(total))-2*phi*(yt(totalminus1)'*yt(total))+((T-1)*(cplusdiff^2))+2*phi*cplusdiff*(ones(167,1)'*yt(totalminus1))+(phi^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtc=(1/0.000001)*(lthetacplussdiff-ltheta);

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

%%%%%%%%%%计算ltheta对phi的导数%%%%%%%%%%%%

phiplusdiff=phi+0.000001;

lthetaphiplusdiff=-0.5*log(2*pi)-0.5*log(sigma2/(1-phiplusdiff^2))-((yt(1)-c/(1-phiplusdiff))^2)/(2*sigma2/(1-phiplusdiff^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigma2)-(1/(2*sigma2))*((yt(total)'*yt(total))-2*c*(ones(167,1)'*yt(total))-2*phiplusdiff*(yt(totalminus1)'*yt(total))+((T-1)*(c^2))+2*phiplusdiff*c*(ones(167,1)'*yt(totalminus1))+(phiplusdiff^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtphi=(1/0.000001)*(lthetaphiplusdiff-ltheta);

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

%%%%%%%%%%计算ltheta对sigma2的导数%%%%%%%%%%%%

sigma2plusdiff=sigma2+0.000001;

lthetasigma2plusdiff=-0.5*log(2*pi)-0.5*log(sigma2plusdiff/(1-phi^2))-((yt(1)-c/(1-phi))^2)/(2*sigma2plusdiff/(1-phi^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigma2plusdiff)-(1/(2*sigma2plusdiff))*((yt(total)'*yt(total))-2*c*(ones(167,1)'*yt(total))-2*phi*(yt(totalminus1)'*yt(total))+((T-1)*(c^2))+2*phi*c*(ones(167,1)'*yt(totalminus1))+(phi^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtsigma2=(1/0.000001)*(lthetasigma2plusdiff-ltheta);

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

% 故在该theta时的梯度为(lthetadiffwrtc,lthetadiffwrtphi,lthetadiffwrtsigma2)'

% 下面求其Hessian

% lthetadiffwrtc再对c求导，表示c发生diff大小变化时，lthetadiffwrtc的变化比上diff

% 用lthetadiffwrtccplusdiff表示c变化diff时的ltheta，则lthetadiffwrtc的变化为

% (lthetadiffwrtccplusdiff-lthetadiffwrtc),则若用lthetadiffwrtcc表示第一行第一个元素

% lthetadiffwrtcc=(1/0.000001)*(lthetadiffwrtccplusdiff-lthetadiffwrtc)

% 注意上一步的运算中，lthetadiffwrtc已经固定了，因此只需要计算lthetadiffwrtccplusdiff

gtheta(1,1)=lthetadiffwrtc;

gtheta(2,1)=lthetadiffwrtphi;

gtheta(3,1)=lthetadiffwrtsigma2;

%%%%%%%%%%%%%%%%%%%%%%%%%%% 第一行 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%计算ltheta对c求导再对c求导%%%%%%%%%%%%%%%%%%%%%%

cplusdiff2=cplusdiff+0.000001;

lthetadiffwrtccplusdiff=-0.5*log(2*pi)-0.5*log(sigma2/(1-phi^2))-((yt(1)-cplusdiff2/(1-phi))^2)/(2*sigma2/(1-phi^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigma2)-(1/(2*sigma2))*((yt(total)'*yt(total))-2*cplusdiff2*(ones(167,1)'*yt(total))-2*phi*(yt(totalminus1)'*yt(total))+((T-1)*(cplusdiff2^2))+2*phi*cplusdiff2*(ones(167,1)'*yt(totalminus1))+(phi^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtcc=(1/0.000001)*(lthetadiffwrtccplusdiff-lthetadiffwrtc);

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

%%%%%%%%%%%%%%%%%%%%%计算ltheta对c求导再对phi求导%%%%%%%%%%%%%%%%%%%%%%

phiplusdiff2=phi+0.000001;

lthetacplussdiffphiplusdiff=-0.5*log(2*pi)-0.5*log(sigma2/(1-phiplusdiff2^2))-((yt(1)-cplusdiff/(1-phiplusdiff2))^2)/(2*sigma2/(1-phiplusdiff2^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigma2)-(1/(2*sigma2))*((yt(total)'*yt(total))-2*cplusdiff*(ones(167,1)'*yt(total))-2*phiplusdiff2*(yt(totalminus1)'*yt(total))+((T-1)*(cplusdiff^2))+2*phiplusdiff2*cplusdiff*(ones(167,1)'*yt(totalminus1))+(phiplusdiff2^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtcphi=(1/0.000001)*(lthetacplussdiffphiplusdiff-lthetadiffwrtc);

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

%%%%%%%%%%%%%%%%%%%%%计算ltheta对c求导再对sigma2求导%%%%%%%%%%%%%%%%%%%%%%

sigmaplusdiff2=sigma2+0.000001;

lthetacplussdiffsigma2plussdiff=-0.5*log(2*pi)-0.5*log(sigmaplusdiff2/(1-phi^2))-((yt(1)-cplusdiff/(1-phi))^2)/(2*sigmaplusdiff2/(1-phi^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigmaplusdiff2)-(1/(2*sigmaplusdiff2))*((yt(total)'*yt(total))-2*cplusdiff*(ones(167,1)'*yt(total))-2*phi*(yt(totalminus1)'*yt(total))+((T-1)*(cplusdiff^2))+2*phi*cplusdiff*(ones(167,1)'*yt(totalminus1))+(phi^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtcsigma2=(1/0.000001)*(lthetacplussdiffsigma2plussdiff-lthetadiffwrtc);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%% 第二行 %%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%计算ltheta对phi求导再对c求导%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% 根据Young定律，lthetadiffwrtphic等于lthetadiffwrtcphi %%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%计算ltheta对phi求导再对phi求导%%%%%%%%%%%%%%%%%%%%%%

phiplusdiff2=phiplusdiff+0.000001;

lthetaphiplusdiffphiplusdiff=-0.5*log(2*pi)-0.5*log(sigma2/(1-phiplusdiff2^2))-((yt(1)-c/(1-phiplusdiff2))^2)/(2*sigma2/(1-phiplusdiff2^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigma2)-(1/(2*sigma2))*((yt(total)'*yt(total))-2*c*(ones(167,1)'*yt(total))-2*phiplusdiff2*(yt(totalminus1)'*yt(total))+((T-1)*(c^2))+2*phiplusdiff2*c*(ones(167,1)'*yt(totalminus1))+(phiplusdiff2^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtphiphi=(1/0.000001)*(lthetaphiplusdiffphiplusdiff-lthetaphiplusdiff);

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

%%%%%%%%%%%%%%%%%%%%%计算ltheta对phi求导再对sigma2求导%%%%%%%%%%%%%%%%%%%%%%

lthetaphiplusdiffsigma2plusdiff=-0.5*log(2*pi)-0.5*log(sigmaplusdiff2/(1-phiplusdiff^2))-((yt(1)-c/(1-phiplusdiff))^2)/(2*sigmaplusdiff2/(1-phiplusdiff^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigmaplusdiff2)-(1/(2*sigmaplusdiff2))*((yt(total)'*yt(total))-2*c*(ones(167,1)'*yt(total))-2*phiplusdiff*(yt(totalminus1)'*yt(total))+((T-1)*(c^2))+2*phiplusdiff*c*(ones(167,1)'*yt(totalminus1))+(phiplusdiff^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtphisigma2=(1/0.000001)*(lthetaphiplusdiffsigma2plusdiff-lthetaphiplusdiff);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%% 第三行 %%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%计算ltheta对sigma2求导再对c求导%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% 根据Young定律，lthetadiffwrtsigma2c等于lthetadiffwrtcsigma2 %%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%计算ltheta对sigma2求导再对phi求导%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% 根据Young定律，lthetadiffwrtsigma2phi等于lthetadiffwrtphisigma2 %%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%计算ltheta对sigma2求导再对sigma2求导%%%%%%%%%%%%%%%%%%%%%%

sigma2plusdiff2=sigma2plusdiff+0.000001;

lthetasigma2plusdiffsigma2plusdiff=-0.5*log(2*pi)-0.5*log(sigma2plusdiff2/(1-phi^2))-((yt(1)-c/(1-phi))^2)/(2*sigma2plusdiff2/(1-phi^2))-((T-1)/2)*log(2*pi)-((T-1)/2)*log(sigma2plusdiff2)-(1/(2*sigma2plusdiff2))*((yt(total)'*yt(total))-2*c*(ones(167,1)'*yt(total))-2*phi*(yt(totalminus1)'*yt(total))+((T-1)*(c^2))+2*phi*c*(ones(167,1)'*yt(totalminus1))+(phi^2)*(yt(totalminus1)'*yt(totalminus1)));

lthetadiffwrtsigma2sigma2=(1/0.000001)*(lthetasigma2plusdiffsigma2plusdiff-lthetadiffwrtcsigma2);

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

% 构造Hessian

hessiantheta(1,1)=lthetadiffwrtcc;

hessiantheta(1,2)=lthetadiffwrtcphi;

hessiantheta(1,3)=lthetadiffwrtcsigma2;

hessiantheta(2,1)=lthetadiffwrtcphi;

hessiantheta(2,2)=lthetadiffwrtphiphi;

hessiantheta(2,3)=lthetadiffwrtphisigma2;

hessiantheta(3,1)=lthetadiffwrtcsigma2;

hessiantheta(3,2)=lthetadiffwrtphisigma2;

hessiantheta(3,3)=lthetadiffwrtsigma2sigma2;

% 设定tolerant

% 更新theta

newtheta=theta+inv(hessiantheta)*gtheta;

newc=newtheta(1,1);

newphi=newtheta(2,1);

newsigma2=newtheta(3,1);

diffc=abs(newc-c);

diffphi=abs(newphi-phi);

diffsigma2=abs(newsigma2-sigma2);

c=newc;phi=newphi;sigma2=newsigma2;