//AR(3)
new;
cls;
library pgraph;  /* present graphs */
pqgwin auto;     /* present plural graphs */
load rgdp[208,1]  = us_rgdp.txt;     //1960q1~2011q4, 208 obs
load e_a[92,1]    = resid_6083.txt;  //1960q1~1983q4, 92  obs
load e_b[112,1]   = resid_8411.txt;  //1984q1~2011q4, 112 obs
lnrgdp  =  ln(rgdp);                       //1960q1~2011q4, 208 obs
dlnrgdp =  lnrgdp[2:208]-lnrgdp[1:208-1];  //1960q2~2011q4, 207 obs
//since AR(3) loses first three data points, we have to start from 1961q1
y  =  lnrgdp[5:208];     //1961q1~2011q4, 205 obs
dy =  dlnrgdp[4:208-1];  //1961q1~2011q4, 205 obs
t=rows(y);               //t=204
R=1|0|0;
/* a.pre 1984 sample */
//y_t=mu_a + phi_1a * y_t-1 + phi_2a * y_t-2 + phi_3a * y_t-3 + e_t , e_t ~ iid (0, sig^2_a)
T_a      =   92;
mu_a     =   0.009036;
phi_1a   =   0.243911;
phi_2a   =   0.113029;
phi_3a   =   -0.030051;
sig2_a   =   0.009392;  // Parameters are from eviews
phi1a    =   1 - phi_1a - phi_2a - phi_3a;
psi1a    =   phi1a;
"psi(1) (pre 1984) is";; psi1a;

y_a  =  y[1:T_a];    //1961q1~1983q4, 92 obs
dy_a =  dy[1:T_a];   //1961q1~1983q4, 92 obs

yz_a   =   y[3:T_a];    //1961q3~1983q4 90 obs
yt_a   =   y[2:T_a-1];  //1961q2~1983q3 90 obs
yl_a   =   y[1:T_a-2];  //1961q1~1983q2 90 obs
dyz_a  =   dy[3:T_a];   //1961q3~1983q4 90 obs
dyt_a  =   dy[2:T_a-1]; //1961q2~1983q3 90 obs
dyl_a  =   dy[1:T_a-2]; //1961q1~1983q2 90 obs

/*
et_a=e_a[2:T_a];        //1961Q2~1983Q4, 91 obs
el_a=e_a[1:T_a-1];      //1961Q1~1983Q3, 91 obs
*/

dy_ztl_a = dyz_a~dyt_a~dyl_a;  // 90 x 3 Matrix (not 3 x 90 ), 1961q3~1983q4, 90 obs
mu_til_a = mu_a|0|0;           //3 x 1 vector
    F_a=
    (phi_1a ~ phi_2a ~ phi_3a)|
    (1      ~ 0      ~ 0     )|
    (0      ~ 1      ~ 0     );
/* B. Post 1984 Sample */
//y_t=mu_b + phi_1b * y_t-1 + phi_2b * y_t-2 + phi_3b * y_t-3 + e_t , e_t ~ iid (0, sig^2_b)
T_b       =   112;
mu_b      =   0.006593;
phi_1b    =   0.378274;
phi_2b    =   0.314645;
phi_3b    =   -0.125415;  
sig2_b    =   0.003083; // Parameters are from Eviews
phi1b   =  1 - phi_1b - phi_2b;
psi1b   =  phi1b;
"psi(1) (post 1984) is";; psi1b;

    F_b = (phi_1b ~ phi_2b ~ phi_3b)|
          (1      ~ 0      ~ 0     )|
          (0      ~ 1      ~ 0     );
         
y_b   =  y[T_a+1:T];   //1984Q1~2011Q4, 112 obs
dy_b  =  dy[T_a+1:T];  //1984Q1~2011Q4, 112 obs
yz_b  = y[T_a+3:T];    //1984q3~2011q4, 110 obs
yt_b  = y[T_a+2:T-1];  //1984Q2~2011Q3, 110 obs
yl_b  = y[T_a+1:T-2];  //1984Q1~2011Q2, 110 obs
dyz_b = dy[T_a+3:T];   //1984q3~2011q4, 110 obs
dyt_b = dy[T_a+2:T-1]; //1984Q2~2011Q3, 110 obs
dyl_b = dy[T_a+1:T-2]; //1984Q1~2011Q2, 110 obs
/*
et_b = e_b[2:T_b];     //1984Q2~2011Q4, 111 obs
el_b = e_b[1:T_b-1];   //1984Q1~2011Q3, 111 obs
*/
dy_ztl_b = dyz_b~dyt_b~dyl_b; //110 x 3 Matrix (not 4 x 110 ), 1984Q3~2011Q4, 110 obs
mu_til_b = mu_b|0|0;//3 x 1 vector
   

@====Impulse Response Analysis : dy_t===============================@
// A. Pre 1984
IR_dy_a=zeros(T_a,1);
    i=1;
    F=eye(3);
   
    do until i>T_a;
   
        FR = F * R;
        IR_dy_a[i,1] = FR[1,1];
        F = F * F_a;
        
    i=i+1;
   
    endo;
title("IR_dy : Pre 1984");
xy(t_a,IR_dy_a); //Graph
// B. Post 1984
IR_dy_b = zeros(T_b,1);
    i =1;
    F=eye(3);
   
    do until i>T_b;
   
        FR=F*R;
        IR_dy_b[i,1]=FR[1,1];
        F=F*F_b;
   
    i=i+1;
   
    endo;
title("IR_dy : Post 1984");
xy(t_b,IR_dy_b); //Graph
@==================================================================@
@====Impulse Response Analysis : y_t===============================@

// A. Pre 1984
IR_y_a=zeros(T_a,1);
    i=1;
   
    do until i>T_a;
   
        IR_y_a[i,1]=sumc(IR_dy_a[1:i]);
   
    i=i+1;
   
    endo;
title("IR_y : Pre 1984");
xy(t_a,IR_y_a); //Graph

// B. Post 1984
IR_y_b=zeros(T_b,1);
    i=1;
   
    do until i>T_b;
   
        IR_y_b[i,1]=sumc(IR_dy_b[1:i]);
   
    i=i+1;
   
    endo;
title("IR_y : Post 1984");
xy(t_b,IR_y_b); //Graph

@============================================================@
@================Prediction : E(dy_t+j|I_t)==================@
xtics(0,100,10,5);
Ytics(-0.03,0.05,0.01,0.001);
// A. Pre 1984
m=10;// maximum prediction horizon
j=1; // Prediction Horizon
prd_dy_a = dyz_a; //1961Q3~1983Q4, 90 obs
   do until j>m;

           prd=zeros(T_a-1,1);
           
            i=1;
            
            do until i>T_a-j-1;
               
                pred=(eye(3)-(F_a^j))*inv(eye(3)-F_a)*mu_til_a+(F_a^j)*(dy_ztl_a[i,.])';//E(dy_til_i+j|I_i)
                prd[i+j,1]=pred[1,1];// The ith element of prd is E(dy_i+j|I_i)
            
            i=i+1;
            
            endo;
            
    prd_dy_a=prd_dy_a~prd;
    j=j+1;     
   
    endo;

title("Prediction E(dy_t+1|I_t) : Pre 1984");
xy(t,prd_dy_a[3:T_a-1,2]~prd_dy_a[3:T_a-1,1]);
title("Prediction E(dy_t+3|I_t) : Pre 1984");
xy(t,prd_dy_a[4:T_a-1,4]~prd_dy_a[4:T_a-1,1]);
title("Prediction E(dy_t+6|I_t) : Pre 1984");
xy(t,prd_dy_a[7:T_a-1,7]~prd_dy_a[7:T_a-1,1]);
title("Prediction E(dy_t+10|I_t) : Pre 1984");
xy(t,prd_dy_a[11:T_a-1,11]~prd_dy_a[11:T_a-1,1]);


红字部分，是ARMA(2，2)的内容命令，怎么样改成AR(3)的命令，xy(t,prd_dy_a[3:T_a-1,2]~prd_dy_a[3:T_a-1,1]);这个表示的是什么意思，Y轴怎么变成矩阵了？？