Gauss 编程应用于学科前沿一：单位根

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Gauss 编程应用于学科前沿一：单位根
单位根理论我就不说了，虽然已经有30余年的历史了，发展也很成熟，但做实证特别是时间序列和长序列面板数据模型时，还需要用到单位根检验。
因此，本帖子从最初的单位根检验开始，到非线性单位根检验，再到面板数据单位根，我们希望通过本帖，基本上掌握单位根检验的各种程序。

在code里，我们有带实际数据的，也有模拟数据，当然也可能只有调用程序，都能用。

我们先从DF检验开始：
1、这个例子是DF检验的模拟，重复10000次，并计算临界值，注意临界值是在原假设下估计的，也就是先生成带单位根的数据AR，即proc gen_ar1(phi0,e_mat); AR程序你也可以自己编，不像这样。

/* Dickey Fuller Distribution */
   new;cls;
   library pgraph;
   tru_phi=1;
   H0_phi=1;
   n_ob =100;
   n_trial =10000;
   t_mat=zeros(n_trial,1);
trial=1;
do until trial > n_trial;
   trial;;"th Iteration";
   e=rndn(n_ob,1);
   y=gen_ar1(tru_phi,e);
   Yt = y[2:n_ob,1];  
    Ylag=y[1:n_ob-1,1];

   tt=rows(Yt);
  X=Ylag~ones(tt,1)~seqa(1,1,tt) ;
   b=inv(X'X)* X'Yt;
   residual=Yt-X*b;
   sigma2=residual'residual/(tt-cols(x));
   vc_b = sigma2 * inv(X'X);   @Variance-Covariance matrix of b@
   se_b =sqrt(vc_b[1,1]);  @Standard error of b@
   t_value=(b[1,1]-H0_phi)/se_b;  @test statistic for H0: phi= H0_phi@
   t_mat[trial,1]=t_value;
trial=trial +1;
endo;
sorted_t = sortc(t_mat,1);
output file =dky_flr.out reset;
    "5% Critical Value is:" sorted_t[500,1];
output off;
{a,b,c}= histp(t_mat,50);
end;
@======================================@
proc gen_ar1(phi0,e_mat);
     local x0_mat,i;
     x0_mat = zeros(n_ob,1);
     x0_mat[1,1]=e_mat[1,1];
     i=2;
     do until i>n_ob;
         x0_mat[i,1]=phi0*x0_mat[i-1,1] + e_mat[i,1];
     i=i+1;
     endo;
retp(x0_mat);
endp;

扫码加我 拉你入群

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

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

分享0 收藏4 回帖

关键词：GAUSS 学科前沿 Aus USS distribution 编程

再看一例，上面的解释变量带常数项和趋势项，
下面我们不考虑这些。
/*========= ============================
Test statistics = sqrt(T)*(phi_hat - H0_phi)
DGP : Y(t) =phi*Y(t-1) + Et,    phi=0
Reg : Y(t) =b1*Y(t-1) +Et
=======================================*/
new;cls;
library pgraph;  @调用图函数库@
   tru_phi=0;
   H0_phi=0;
   n_ob =100;
   n_trial =1000;
   count=0;
   t_mat=zeros(n_trial,1);

ii=1;
do until ii> n_trial;
     trial;;"th Iteration";
@============= DGP =============@
     ymat = zeros(n_ob,1);
     ymat[1]=rndn(1,1)*sqrt((1/(1-tru_phi^2)));
     i=2;
     do until i>n_ob;
         ymat[i,1]=tru_phi*ymat[i-1,1] + rndn(1,1);
     i=i+1;
     endo;
@===========回归==============@
   Yt = ymat[2:n_ob,1];  
   Ylag=ymat[1:n_ob-1,1];
   X=Ylag;

   b=inv(X'X)* X'Yt;

@==========Test Statistic==========@
   static=sqrt(n_ob)*(b-H0_phi);                @test statistic for H0: phi= H0_phi@
   t_mat[ii,1]=static;
@====Counting the # of rejction H0===@
   if abs(static) > 1.96;  @Critical value at a 5% significance level@
       count=count+1;           
   endif;
ii=ii +1;
endo;

output file = dky_flr.prn reset; t_mat; output off;
{a,b,c}= histp(t_mat,50);   @直方图@
output file = dky_flr.out reset;
"Number of rejections from conventional t-test ";   count;
output off;
end;

如果回归方程带常数项 : Y(t) =b0+b1*Y(t-1) +Et

则把上面程序中的x换为：X=ones(n_ob-1,1)~Ylag;就可以了

上面两例如果大家注意观看图型，我们会发现基本上是正态分布，因为我们设置了tru_phi=0;
下面我使用不同的统计量，使用tru_phi=1; 该统计量是偏态，也就是说当有单位根了，传统的渐近分布就不能用了。
读者可以看陆懋祖的高等计量经济学，有关单位根的渐近分布推导。

/*========= ============================
检验统计量不是根据t统计量计算的，现在是： T*(phi_hat - H0_phi)
DGP : Y(t) =phi*Y(t-1) + Et,  phi=1
Reg : Y(t) =b1*Y(t-1) +Et
=======================================*/
new;cls;
library pgraph;
   tru_phi=1;
   H0_phi=1;
   n_ob =100;
   n_trial =1000;
   count=0;
   t_mat=zeros(n_trial,1);

trial=1;
do until trial > n_trial;
   trial;;"th Iteration";
@============= DGP =============@
     ymat = zeros(n_ob,1);
     i=2;
     do until i>n_ob;
         ymat[i,1]=tru_phi*ymat[i-1,1] + rndn(1,1);
     i=i+1;
     endo;
@=============Reg==============@
   Yt = ymat[2:n_ob,1];  Ylag=ymat[1:n_ob-1,1];
   X=Ylag;
   b=inv(X'X)* X'Yt;

@==========Test Statistic==========@
   static=(n_ob)*(b-H0_phi);  @test statistic for H0: phi= H0_phi@
   t_mat[trial,1]=static;
@====Counting the # of rejction H0===@
   if static < -8;    @Testing the hypothesis at a 5% significance
                                   level = -8 for case 1 @
        count=count+1;
   endif;

trial=trial +1;
endo;

{a,b,c}= histp(t_mat,50);
output file = dky_flr.out reset;

"Number of rejections from conventional t-test ";   count;
output off;
end;

读者最好根据例一中的临界值程序得到相关的临界值，然后在对后面例子中的
1.96进行替代，这样就更精确了，比如你的数据是40个，这个时候可以模拟得到40个数据的临界值，虽然相差不大，但很有用。
   if abs(static) > 1.96;        @Testing the hypothesis at a 5% significance
                                   level@
        count=count+1;
   endif;

明天继续，希望用1个月的时间把单位根相关问题都整理处理。

好东西，打包下就好了，希望gauss板块活跃起来

df.txt (1.77 KB)

2012-3-16 18:21:31 上传

数据

继续昨天的工作。
本例使用实际数据做DF检验，数据在附件
new; cls;
load data[50,2] = df.txt;
y=data[.,1];     @数据@
n=rows(data);
dy=y[2:n]-y[1:n-1];  @y差分@
y_1=y[1:n-1];           @原检验方程为：yt=rho*yt_1+et，检验H0: rho=1 ；也可以变为：dy=delta*yt_1+et ，检验H0: delta=0，y_1为解释变量@
call dftau(dy,y_1);
proc dftau(dy,x);
local n1,k,b,e,s2hat,varb,se,tau;
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1];
   print "tau=" tau;; print/rz "( n=" n1 ")";
   retp(tau);
endp;

模拟DF检验的分位数
new; cls;
n=50;                   @修改样本量，可以得到各种样本下的分位数@
pp={0.01,0.05,0.10};  @3个百分位数@
times=10000;             @循环次数@
tau=zeros(times,1);     
   j=1;
   do while j<=times;
      e=rndn(n+1,1);
      y=zeros(n+1,1);
      y0=0;
      y[1]=y0+e[1];
      i=2;
      do while i<=n+1;
         y=y[i-1]+e;
         i=i+1;
      endo;
      n1=rows(y);
      dy=y[2:n1]-y[1:n1-1];
      y_1=y[1:n1-1];
      tau[j]=dftau(dy,y_1);
      j=j+1;
   endo;
   tau=sortc(tau,1);
   stat=tau[round(pp*times)];         @round表示取整@
print stat;

proc dftau(dy,x);
   local n1,k,b,e,s2hat,varb,se,tau;
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1];
   retp(tau);
endp;

DF检验有三种形式，对于带常数项和趋势项的见后两个：
I:   dy=delta+y_1+et
II:  dy=delta*y_1+alpha+et
III: dy=delta*y_1+alpha+beta*t+et
对于第二种形式的DF检验，还可以用F统计量。
H0: delta=0     H1: delta<0
限制alpha=0，和delta=0，这样有约束R*beta=q；
   R={1 0,
      0 1}
   q={0,
      0}.
非约束模型是dy=delta*y_1+alpha+et；约束模型是 alpha=0，和delta=0。则受约束的系数估计为
br=b+inv(x'x)*R'inv(R*inv(x'x)*R')*(q-R*b)
计算F统计量F=((RSSr-RSSu)/J)/(RSSu/(n-k))；RSSr约束下的残差平方和，
RSSu非约束下的。

new;cls;
load data[50,1]=df.txt;
y=data[.,1];
n=rows(y);
dy=y[2:n]-y[1:n-1];
y_1=y[1:n-1];
call dfF2(dy,y_1);
proc dfF2(dy,x);
   local n1,k,b,rssu,rssr,R,q,br,fstat;
   n1=rows(x); k=2;
   x=ones(n1,1)~x;
   b=inv(x'x)*x'dy;
   rssu=(dy-x*b)'(dy-x*b);
   R={1 0,
      0 1};
   q={0,
      0};
   br=b+inv(x'x)*R'inv(R*inv(x'x)*R')*(q-R*b);
   rssr=(dy-x*br)'(dy-x*br);
   fstat=((rssr-rssu)/2)/(rssu/(n1-k));
   print "F=" fstat;; print/rz "( n=" n1 ")";
   retp(fstat);
endp;