@===========================================================================@
new;
library optmum, pgraph;
format /m1 /rd 8,4;
load data[191,3]=s&w.prn; @1995.1 -- 2010.11@
@ month, yiban, jiagong@
lndata=ln(data);
month=data[2:191,1];
y1=(ln(data[2:191,2]) - ln(data[1:190,2]))*100;
y2=(ln(data[2:191,3]) - ln(data[1:190,3]))*100;
yy0=(y1~y2);
@================================Option 2: Based on======================
Demeaned data
========================================================================@
yy_m=meanc(yy0);
yy=yy0'-yy_m; @2xT@
prmtr_in={0.2 0.1 0.1 0.1
0.2 0.1 0.1 0.2 0.1 0.1
0.5 0.5
0.5 0 0 0.5 0 0 };
@=====================================================================@
t=cols(yy);
@========Maximum Likelihood Estimation=========================================@
prmtr_in=prmtr_in';
{xout,fout,gout,cout}=optmum(&lik_fcn,prmtr_in);
@========End of Main Program===================================================@
proc lik_fcn(prmtr1);
local sig_e,sig,F,H,Q,R,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT,vt,ft,
valu,loglik,j_iter,vecP_LL,prmtr,phi1,phi2,phi3,phi4,psi11,psi12,psi13,
psi21,psi22,psi23,sig1,sig2,val,gamma11,
gamma12,gamma13,gamma21,gamma22,gamma23;
loglik=zeros(t,1);
prmtr=trans(prmtr1);
locate 20,1; prmtr';
phi1=prmtr[1,1];
phi2=prmtr[2,1];
phi3=prmtr[3,1];
phi4=prmtr[4,1];
psi11=prmtr[5,1];
psi12=prmtr[6,1];
psi13=prmtr[7,1];
psi21=prmtr[8,1];
psi22=prmtr[9,1];
psi23=prmtr[10,1];
sig1=prmtr[11,1];
sig2=prmtr[12,1];
gamma11=prmtr[13,1];
gamma12=prmtr[14,1];
gamma13=prmtr[15,1];
gamma21=prmtr[16,1];
gamma22=prmtr[17,1];
gamma23=prmtr[18,1];
@=====================================================================@
F=(phi1~ phi2~ phi3~ phi4~ 0~ 0~ 0~ 0~ 0~ 0)|
( 1~ 0~ 0~ 0~ 0~ 0~ 0~ 0~ 0~ 0)|
( 0~ 1~ 0~ 0~ 0~ 0~ 0~ 0~ 0~ 0)|
( 0~ 0~ 1~ 0~ 0~ 0~ 0~ 0~ 0~ 0)|
( 0~ 0~ 0~ 0~ psi11~ psi12~ psi13~ 0~ 0~ 0)|
( 0~ 0~ 0~ 0~ 1~ 0~ 0~ 0~ 0~ 0)|
( 0~ 0~ 0~ 0~ 0~ 1~ 0~ 0~ 0~ 0)|
( 0~ 0~ 0~ 0~ 0~ 0~ 0~ psi21~ psi22~ psi23)|
( 0~ 0~ 0~ 0~ 0~ 0~ 0~ 1~ 0~ 0)|
( 0~ 0~ 0~ 0~ 0~ 0~ 0~ 0~ 1~ 0);
H=(gamma11~ gamma12~ gamma13~ 0~ 1~ 0~ 0~ 0~ 0~ 0)|
(gamma21~ gamma22~ gamma23~ 0~ 0~ 0~ 0~ 1~ 0~ 0);
sig_e=1;
Q=(sig_e~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ sig1~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ sig2~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0);
R=0;
@================================================================@
BETA_LL=0|0|0|0|0|0|0|0|0|0;
vecP_LL=inv(eye(10^2) - F.*.F)*vec(Q); @..Vectorization..@
P_LL=vecP_LL[1:10,1]~vecP_LL[11:20,1]~vecP_LL[21:30,1]~vecP_LL[31:40,1]
~vecP_LL[41:50,1]~vecP_LL[51:60,1]~vecP_LL[61:70,1]~vecP_LL[71:80,1]
~vecP_LL[81:90,1]~vecP_LL[91:100,1];
@..initial state vector and its covariance ..@
j_iter=1;
do until j_iter > t;
BETA_TL = F*BETA_LL;
P_TL = F*P_LL*F' + Q;
vt = yy[.,j_iter] - H*BETA_TL;
ft = H*P_TL*H' + R;
BETA_TT = BETA_TL + P_TL*H'*inv(ft)*vt;
P_TT = P_TL - P_TL*H'*inv(ft)*H*P_TL;
BETA_LL=BETA_TT;
P_LL=P_TT;
val=v_prob(vt,ft); @likelihood@
loglik[j_iter,1] = -1*ln(val); @ - log likelihood@
j_iter=j_iter+1;
endo;
retp(sumc(loglik));
endp;
@==============================================================================@
@==============================================================================@
proc (2)= filter(prmtr1);
local sig_e,sig,F,H,Q,R,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT,vt,ft,
valu,loglik,j_iter,vecP_LL,prmtr,phi1,phi2,phi3,phi4,psi11,psi12,psi13,
psi21,psi22,psi23,sig1,sig2,val,gamma11,
gamma12,gamma13,gamma21,gamma22,gamma23,out_mat,c_bar,ss_kg,k,alpfa;
prmtr=trans(prmtr1);
phi1=prmtr[1,1];
phi2=prmtr[2,1];
phi3=prmtr[3,1];
phi4=prmtr[4,1];
psi11=prmtr[5,1];
psi12=prmtr[6,1];
psi13=prmtr[7,1];
psi21=prmtr[8,1];
psi22=prmtr[9,1];
psi23=prmtr[10,1];
sig1=prmtr[11,1];
sig2=prmtr[12,1];
gamma11=prmtr[13,1];
gamma12=prmtr[14,1];
gamma13=prmtr[15,1];
gamma21=prmtr[16,1];
gamma22=prmtr[17,1];
gamma23=prmtr[18,1];
@=====================================================================@
F=(phi1~ phi2~ phi3~ phi4~ 0~ 0~ 0~ 0~ 0~ 0)|
( 1~ 0~ 0~ 0~ 0~ 0~ 0~ 0~ 0~ 0)|
( 0~ 1~ 0~ 0~ 0~ 0~ 0~ 0~ 0~ 0)|
( 0~ 0~ 1~ 0~ 0~ 0~ 0~ 0~ 0~ 0)|
( 0~ 0~ 0~ 0~ psi11~ psi12~ psi13~ 0~ 0~ 0)|
( 0~ 0~ 0~ 0~ 1~ 0~ 0~ 0~ 0~ 0)|
( 0~ 0~ 0~ 0~ 0~ 1~ 0~ 0~ 0~ 0)|
( 0~ 0~ 0~ 0~ 0~ 0~ 0~ psi21~ psi22~ psi23)|
( 0~ 0~ 0~ 0~ 0~ 0~ 0~ 1~ 0~ 0)|
( 0~ 0~ 0~ 0~ 0~ 0~ 0~ 0~ 1~ 0);
H=(gamma11~ gamma12~ gamma13~ 0~ 1~ 0~ 0~ 0~ 0~ 0)|
(gamma21~ gamma22~ gamma23~ 0~ 0~ 0~ 0~ 1~ 0~ 0);
sig_e=1;
Q=(sig_e~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ sig1~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ sig2~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0)|
( 0~ 0~0~0~ 0~0~0~ 0~ 0~0);
R=0;
@================================================================@
BETA_LL=0|0|0|0|0|0|0|0|0|0;
vecP_LL=inv(eye(10^2) - F.*.F)*vec(Q); @..Vectorization..@
P_LL=vecP_LL[1:10,1]~vecP_LL[11:20,1]~vecP_LL[21:30,1]~vecP_LL[31:40,1]
~vecP_LL[41:50,1]~vecP_LL[51:60,1]~vecP_LL[61:70,1]~vecP_LL[71:80,1]
~vecP_LL[81:90,1]~vecP_LL[91:100,1];
out_mat=zeros(t,3);
j_iter=1;
do until j_iter > t;
BETA_TL = F*BETA_LL;
P_TL = F*P_LL*F' + Q;
vt = yy[.,j_iter] - H*BETA_TL;
ft = H*P_TL*H' + R;
BETA_TT = BETA_TL + P_TL*H'*inv(ft)*vt;
P_TT = P_TL - P_TL*H'*inv(ft)*H*P_TL;
BETA_LL=BETA_TT;
P_LL=P_TT;
out_mat[j_iter,.]=vt'~beta_tt[1,1];
j_iter=j_iter+1;
endo;
ss_kg=P_TL*H'*inv(ft); @..Steady-State Kalman Gain..@
k=(eye(10)-ss_kg*h)*f;
alpfa=inv(eye(10) - k)*ss_kg*yy_m;
c_bar=alpfa[1,1];
retp(out_mat, c_bar); @Returns out_mat and the mean of the common component@
endp;
@==============================================================================@
proc v_prob(err,var); @..Calculates Pr[Yt/St, Yt-1]..@
local valu;
valu = (1/sqrt(det(var)))*exp(-0.5*err'inv(var)*err);
retp(valu);
endp;
@==============================================================================@
proc trans(c0); @..Constraining Values of Reg. Coff..@
local c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,cc;
cc=c0;
c1=c0[1,1]/(1+abs(c0[1,1]));
c2=c0[2,1]/(1+abs(c0[2,1]));
c3=c0[3,1]/(1+abs(c0[3,1]));
c4=c0[4,1]/(1+abs(c0[4,1]));
cc[1,1]=c1+c2+c3+c4;
cc[2,1]=-1*c1*c2-1*c1*c3-1*c1*c4-1*c2*c3-1*c2*c4-1*c3*c4;
cc[3,1]=c1*c2*c3+c1*c2*c4+c1*c3*c4+c2*c3*c4;
cc[4,1]=-1*c1*c2*c3*c4;
c5=c0[5,1]/(1+abs(c0[5,1]));
c6=c0[6,1]/(1+abs(c0[6,1]));
c7=c0[7,1]/(1+abs(c0[7,1]));
cc[5,1]=c5+c6+c7;
cc[6,1]=-1*c5*c6-c5*c7-c6*c7;
cc[7,1]=c5*c6*c7;
c8=c0[8,1]/(1+abs(c0[8,1]));
c9=c0[9,1]/(1+abs(c0[9,1]));
c10=c0[10,1]/(1+abs(c0[10,1]));
cc[8,1]=c8+c9+c10;
cc[9,1]=-1*c8*c9-c8*c10-c9*c10;
cc[10,1]=c8*c9*c10;
cc[11:12]=c0[11:12]^2;
cc[13:18]=c0[13:18]/10;
retp(cc);
endp;
错误提示如下:
Stack trace:
lik_fcn called from c:\gauss8.0\src\optutil.src, line 151
_optmum called from c:\gauss8.0\src\optmum.src, line 252
optmum called from C:\Documents and Settings\admin\桌面\s&w1.opt, line 51
请教这是什么意思啊?