全部版块 我的主页
论坛 计量经济学与统计论坛 五区 计量经济学与统计软件 Gauss专版
4733 11
2011-03-30
@===========================================================================@
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
请教这是什么意思啊?
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

全部回复
2011-3-30 20:33:48
是不是保存(或调用)路径的问题啊?----弱弱的问问,这是什么语言的程序啊?感觉是SAS,但觉得又不是.....别拍我啊。。。
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-31 09:37:08
Gauss专版,当然是Gauss程序了。闹心啊,改了一个礼拜的程序就卡在这了
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-31 09:37:26
毕业论文等着用呢
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-31 10:39:30
你这是做什么的程序?我运行错误提示如下:

   C:\010\saa.prg(92) : error G0047 : Rows don't match
P_LL=vecP_LL[1:10,1]~vecP_LL[10:20,1]~vecP_LL[20:30,1]~vecP_LL[30:40,1]
        ~vecP_LL[40:50,1]~vecP_LL[50:60,1]~vecP_LL[60:70,1]~vecP_LL[70:80,1]
        ~vecP_LL[80:90,1]~vecP_LL[90:100,1];
                           @..initial state vector and its covariance ..@
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2011-3-31 10:44:32
0.4394  -0.0702   0.0049  -0.0001   0.3485  -0.0386   0.0014   0.3485  -0.0386   0.0014   0.2500   0.2500   0.0500   0.0000   0.0000   0.0500   0.0000   0.0000
以上是程序运行的结果,但是提示了错误,我估计这结果应该有问题
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

点击查看更多内容…
相关推荐
栏目导航
热门文章
推荐文章

说点什么

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群