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

对于第三个模型：dy=alpha+beta*t+delta*y_1+et
约束
alpha=0, beta=0 and delta=0
   R={1 0 0,
      0 1 0,
      0 0 1}
   q={0,
      0,
      0} .


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 dfF3(dy,y_1);
proc dfF3(dy,x);
   local n1,k,b,rssu,rssr,R,q,br,fstat;
   n1=rows(x); k=3;
   x=ones(n1,1)~seqa(1,1,n1)~x;   @解释变量@
   b=inv(x'x)*x'dy;
   rssu=(dy-x*b)'(dy-x*b);
   R={1 0 0,
      0 1 0,
      0 0 1};
   q={0,
      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)/3)/(rssu/(n1-k));
   print "F=" fstat;; print/rz "( n=" n1 ")";
   retp(fstat);
endp;

顶。。。。。

这个不错，加油～

下面这例是ADF检验，程序里有两种选择滞后值的方法。

new;cls;
maxk = 8;  @数据选择滞后值的设置的最大滞后@
n = 100;
crit_t = 1.645;
p = 1;       @ 1 = with trend, 0 = with constant only, -1 = no constant or trend @
select = 0;  @ 1 = fixed lag, 0 = data dependent lag selection @
y = cumsumc(rndn(n,1));  @生成y序列@
@ or, Read the dat file here. load data[n,1] = data.txt;  @
/* Using a Fixed lag */
if select == 1;
   lag = 4;
   {beta,tval, alpha, tstat} = adf(y,p,lag);   @调用ADF估计程序@
   format /rd/m1 12,4;
   "Using a Fixed lag"
   "Fixed lag    = " lag;
   "T(a-1)       = " alpha;
   "ADF stat     = " tstat;
   "Estimated coefficients :
     (x(t-1), dx(t-1), dx(t-2), ... dx(t-k), const, trend) ";
   beta';
   "t-stat";
   tval';
   goto next;
endif;
/* Using data-driven lag selection procedure */
statkeep = zeros(maxk+1,1);
tvalkeep = zeros(maxk,1);
{beta,tval, alpha, tstat} = adf(y,p,0);
statkeep[1] = tval[1];
j=1;
do while j <= maxk;    @构造数据选择滞后值的一个循环@
   {beta,tval, alpha, tstat} = adf(y,p,j);
   statkeep[j+1] = tval[1];   /* statkeep  for k = 0, 1, 2, ..., maxk */
   tvalkeep[j] = tval[j+1];   /* tvalkeep  for k = 1, 2, ..., maxk */
j = j + 1;
endo;
@ selection using t-test，这个可以看作者的文章，通过abs(tvalkeep[j]) > crit_t来选择@
/* select the lag where the coeff is sig. from maxk to 0. if all coeff are insig., kk = 0  */
j = maxk;
do while j >= 0;
        if (abs(tvalkeep[j]) > crit_t) or (j == 0);
            kk = j;
            j = 0;
        endif;
j = j - 1;
endo;
@ kk = selected lag @
{beta, tval, alpha, tstat} = adf(y,p,kk);
   format /rd/m1 12,4;
   "Using Data dependent selection of the lag";
   "Selected lag = " kk;
   "T(a-1)       = " alpha;
   "ADF stat     = " tstat;
   "Estimated coefficients :
     (x(t-1), dx(t-1), dx(t-2), ... dx(t-k), const, trend) ";
   beta';
   "t-stat";
   tval';
next:
end;
/*** ADF
** Format:  {beta, tval, alpha, tstat} = adf(x,p,l);
** Input:   x     -- time series variable
**          p     -- order of the time-polynomial to include in the
**                   ADF regression.  Set p = -1 for no deterministic
**                   part.
**          l     -- number of lagged changes of x to include in the
**                   fitted regression.
** Output:  alpha       --  estimate of the autoregressive parameter;
**          tstat       --  ADF t-statistic
**          beta        --  coeff. of (x(t-1), dx(t-1), dx(t-2), ... dx(t-k),
                                   constant, time trend)
**          tval        --  t-values
*/
proc (4) = adf(x,p,l) ;
local b,k,z,res,so,var_cov,xx, tval;
local timep,t,m,xmat,nobs,dep,ch,f;
     if (p < -1);
        "Error: p cannot be set < -1";
         retp(0,0,zeros(6,1));
     endif ;
     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;
     nobs    = rows(x);
     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;
     dep     = trimr(x,1,0);
     ch      = diff(x,1);
     k       = 0     ;
     z       = trimr(lagn(x,1),1,0) ;
     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);
     if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ;
     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     b[1,1] = b[1,1] - 1;
     tval    = b ./ sqrt(diag(var_cov));
retp(b,tval,b[1,1],(b[1,1])/sqrt(var_cov[1,1])) ;
endp ;
proc lagn(x,n);
    local y;
    y = shiftr(x', n, (miss(0, 0))');
    retp(y');
endp;
proc diff(x,k) ;
     if ( k == 0) ;
        retp(x) ;
     endif ;
     retp(trimr(x,k,0)-trimr(lagn(x,k),k,0)) ;
endp ;

proc ptrend(p,nobs) ;
local data, t, u, m, timep, it, iit, xmat, invx, beta, resid ;
     u        = ones(nobs,1) ;
     if p > 0 ;
       timep    = zeros(nobs,p) ;
       t        = seqa(1,1,nobs);
       m        = 1 ;
       do while m <= p ;
          timep[.,m] = t^m ;
          m = m + 1 ;
       endo ;
       xmat     = u~timep ;
     else ;
       xmat = u ;
     endif ;
retp(xmat) ;
endp ;

下面是ADF检验统计量的百分位数程序，有了百分数就可以进行检验水平和检验势的模拟了。
可以用百分位数作为临界值进行检验。


new;cls;
/*  Options to choose */
q={0.01,0.025,0.05,0.1,0.9,0.95,0.975,0.99};/* 分位数*/
T=1000;
n=100;/* 100 observations + 50 presample values */
seed1=12456456; /* fix the seed */
maxk = 8;
     @ max # lags allowed in the data-driven procedure @
crit_t = 1.645;
     @ critical value used for significance of lagged terms @
p = 1;
     @ 1 = with trend, 0 = with constant only, -1 = no constant or trend @
select = 0;
     @ 1 = fixed lag, 0 = data dependent lag selection @
y0=0;
tau=zeros(T,1);
ptau=zeros(8,1);
i=1;
do while i<=T;
   y1=zeros(n,1);
   u=rndns(n,1,seed1);
   y1[1,1]=y0+u[1,1];
   j=2;
   do while j<=n;
      y1[j,1]=y1[j-1,1]+u[j,1];
      j=j+1;
   endo;
   y= y1[51:n,.];
   /* Using a Fixed lag */
   if select == 1;
      lag = 4;
      {beta,tval, alpha, tstat} = adf(y,p,lag);
      tau[i,.]=tval[1:1,.];
      format /rd/m1 12,4;
      i=i+1;
      goto next;
    endif;
/* Using data-driven lag selection procedure
   'general to specific' search
    see Ng and Perron (1995, JASA) or Lee (1997, WP) */
@ maxk # regressions @
     statkeep = zeros(maxk+1,1);
     tvalkeep = zeros(maxk,1);
     {beta,tval, alpha, tstat} = adf(y,p,0);
     statkeep[1] = tval[1];
     j=1;
     do while j <= maxk;
        {beta,tval, alpha, tstat} = adf(y,p,j);
        statkeep[j+1] = tval[1];
        tvalkeep[j] = tval[j+1];
        j = j + 1;
     endo;
     /* statkeep  for k = 0, 1, 2, ..., maxk */
     /* tvalkeep  for k = 1, 2, ..., maxk */
@ selection using t-test @
     /* select the lag where the coeff is sig. from maxk to 0.
        if all coeff are insig., kk = 0  */
     j = maxk;
     do while j >= 0;
        if (abs(tvalkeep[j]) > crit_t) or (j == 0);
            kk = j;
            j = 0;
        endif;
        j = j - 1;
     endo;
   @ kk = selected lag @
   {beta, tval, alpha, tstat} = adf(y,p,kk);
    tau[i,.]=tval[1:1,.];
   format /rd/m1 12,4;
   i=i+1;
next:
endo;
pt=tau;
ptau=quantile(tau,q);
print ptau;

在把上例中插入调用程序以及其后面的语句

上面是DF和ADF检验，从下面开始，我们对其他单位根检验进行说明，必要的时候我会把文章附上去的。
一、KPSS检验。
kpss检验特点是非参数检验，需要核函数，而且核函数以及窗宽的选择非常关键，本例选择了Bartlett 。


new; cls;
let data=
              80         10236000
              78         10351000
              71          9371000
              73          9126000
              70          8230000
              71         12111000
              72          5364000
              76          6222000
              76          3952000
              78          3859000
              80          4798000
              80          5133000
              77          4240000
              69         13687000
              69          5248000
              68          4058000
              67          5253000
              66          8963000
              67          6457000
              70          8882000
              69          4134000
              68          3736000
              69          3887000
              69          2095000
              68          1980000
              63          7433000
              63          5487000
              63          3996000
              58         10858000
              60         12095000
              57          7251000
              53         17632000
              55         42481000
              65         16051000
              64          9145000
              62         11264000
              64          4799000
              64          6888000
              66          7093000
              69          6598000
              70          6280000
              69          5372000
              65          7381000
              64          2792000
              63          4258000
              60          3074000
              62          4427000
              60          3767000
              61          3041000
              60         10240000
;
data=reshape(data,rows(data)/2,2);
y=data[.,1];     @数据可以改成自己的@

call kpss(y,4);  @kpss（）中后面的数字是p -- specified max order of autocorrelation@

/***************************************************************************
proc            : KPSS
Input           : y -- (T x 1) vector of obs on variable of interest
                : p -- specified max order of autocorrelation
Output          : Prints KPSS statistics
Notes           : KPSS(trend) based on y(t)=a+bt+e(t) -- trend stationarity
                  KPSS(level) based on y(t)=a+e(t)    -- level stationarity
***************************************************************************/
proc(0)=kpss(y,l);
local T,trend,X,b,e,g0,p,g,i,w,s2,S,ei,etathat,etamhat,lag,eta,oldnv;
    /* Level-stationary test */
    T=rows(y);              @ # of obs @
    X=ones(T,1);            @ RHS var @
    b=y/X;                  @ OLS estimates @
    e=y-X*b;                @ OLS resids @
    g0=e'e/T;               @ gamma_0=(1/T)sum[e(t)^2] @
    etamhat=zeros(l+1,1);
    p=1;
    do until p>l;
        g=zeros(p,1);       @ autocovs @
        i=1;
        do until i>p;
            g[i,1]=(e[1+i:T,1]'e[1:T-i,1])/T;         @ gamma_i=(1/T)sum[e(t)e(t-i)] @
            i=i+1;
        endo;
        w=zeros(p,1);
        i=1;
        do until i>p;
            w[i,1]=(p+1-i)/(p+1);   @ Bartlett window weight @
            i=i+1;
        endo;
        s2=g0+2*w'g;                @ consistent error variance estimate @
        S=zeros(T,1);               @ resid partial sum process @
        i=1;
        do until i>T;
            ei=e[1:i,.];
            S[i,.]=sumc(ei);        @ S(i)=sum[e(i)] @
            i=i+1;
        endo;
        etamhat[1,.]=(1/(g0*T^2))*S'S;      @ KPSS eqn (13), l=0 @
        etamhat[p+1,.]=(1/(s2*T^2))*S'S;    @ KPSS eqn (13) @
        p=p+1;
    endo;
    /* Trend-stationary test */
    trend=seqa(1,1,T);      @ linear time trend @
    X=ones(T,1)~trend;      @ regressor matrix @
    b=y/X;                  @ OLS estimates @
    e=y-X*b;                @ OLS resids @
    g0=e'e/T;               @ gamma_0=(1/T)sum[e(t)^2] @
    etathat=zeros(l+1,1);
    p=1;
    do until p>l;
        g=zeros(p,1);       @ autocovs @
        i=1;
        do until i>p;
            g[i,1]=(e[1+i:T,1]'e[1:T-i,1])/T;@ gamma_i=(1/T)sum[e(t)e(t-i)] @
            i=i+1;
        endo;
        w=zeros(p,1);
        i=1;
        do until i>p;
            w[i,1]=(p+1-i)/(p+1);   @ Bartlett window weight @
            i=i+1;
        endo;
        s2=g0+2*w'g;                @ consistent error variance estimate @
        S=zeros(T,1);               @ resid partial sum process @
        i=1;
        do until i>T;
            ei=e[1:i,.];
            S[i,.]=sumc(ei);        @ S(i)=sum[e(i)] @
            i=i+1;
        endo;
        etathat[1,.]=(1/(g0*T^2))*S'S;      @ KPSS eqn (17), l=0 @
        etathat[p+1,.]=(1/(s2*T^2))*S'S;    @ KPSS eqn (17) @
        p=p+1;
    endo;
    format 8,4;
    "One-sided test of H0: Stationary vs. H1: Unit root";?;
    "Approximate asymptotic critical values";
    "--------------------------------------------------";
    "              10%     5%      1%";
    "--------------------------------------------------";
    "Level        0.347   0.463   0.739";
    "Trend        0.119   0.146   0.216";
    "--------------------------------------------------";?;
    lag=seqa(0,1,l+1);
    eta=etamhat~etathat;
    oldnv=__fmtnv;
    __fmtnv="*.*lg"~8~3;
    "KPSS stats";
    "--------------------------------------------------";
    "      Lag     Level    Trend";
    "--------------------------------------------------";
    i=1;
    do until i>l+1;
        call printfmt(lag[i,.],1);"    " eta[i,.];
        i=i+1;
    endo;
    "--------------------------------------------------";
    "Obs";;call printfmt(T,1);
    __fmtnv=oldnv;
endp;

对于较高的自回归参数，标准的kpss存在过度的检验水平。我们需要利用异方差与自回归一致协方差估计量来估计自回归过程的长期方差，在有限样本下，窗宽选择过大，长期方差过度估计，检验统计量变得非常小，从而检验功效变得很小；选择过小和较大的自回归参数，长期方差会低估，检验统计量变得非常大，有较大的检验水平。Bart Hobijn（1998）利用Newey and West，1994年的方法，即依数据自动选择窗宽的方法，调整了标准kpss的不足。
有兴趣的 kpsstest.pdf (309.62 KB)

2012-4-7 03:22:30 上传

读者可以读读这篇文章。

下面的gauss 程序，考虑了两种核函数，也是各种非参估计中非常重要的两种：Quadratic spectral kernel、Bartlett kernel ，特别是两种核函数在gauss中的设置，而且这些都可以应用到其他的情况，特别是有非参估计的问题。程序有点长。

new; cls;
let data=
              80         10236000
              78         10351000
              71          9371000
              73          9126000
              70          8230000
              71         12111000
              72          5364000
              76          6222000
              76          3952000
              78          3859000
              80          4798000
              80          5133000
              77          4240000
              69         13687000
              69          5248000
              68          4058000
              67          5253000
              66          8963000
              67          6457000
              70          8882000
              69          4134000
              68          3736000
              69          3887000
              69          2095000
              68          1980000
              63          7433000
              63          5487000
              63          3996000
              58         10858000
              60         12095000
              57          7251000
              53         17632000
              55         42481000
              65         16051000
              64          9145000
              62         11264000
              64          4799000
              64          6888000
              66          7093000
              69          6598000
              70          6280000
              69          5372000
              65          7381000
              64          2792000
              63          4258000
              60          3074000
              62          4427000
              60          3767000
              61          3041000
              60         10240000
;
data=reshape(data,rows(data)/2,2);
y=data[.,1];
kerntype = 1;
lagtype = {1,2};
{x0} = wtest( y ,1 ,kerntype ,lagtype );
print x0;
/*** KPSSTest.SRC:
** Wtest            Generalized KPSS test.
** Bartlett         Bartlett spectral window
** QS               Quadratic Spectral window
** Lagsel           Newey and West (1994) automatic lag selection procedure.
** NB               Bandwidth choice for Lagsel using Bartlett.
** NQS              Bandwidth choice for Lagsel using Quadratic Spectral.
** MB               Bandwidth choice for Bartlett without lagsel.
** MQS              Bandwidth choice for Quadratic Spectral without lagsel.
** Autocov          Calculates autocovariances.
*/
/*** Procedure wtest:
** Format:
** w = Wtest( y , h , kerntype , lagtype );
** Input:
** y        (Tx1)-vector, containing time series of which stationarity is tested.
** h        scalar, null hypothesis:
**              = 0     zero mean stationarity
**              = 1     level stationarity
**              = 2     trend stationarity
** kerntype scalar, kernel used:
**              = 0     Bartlett kernel
**              = 1     Quadratic Spectral kernel
** lagtype  (2x1)-matrix: lag selection procedure.
**              lagtype[1]      = 0      no automatic lag selection (m=n)
**                              = 1      automatic lag selection
**              lagtype[2]      scalar/integer, bandwidth. (n)
***/
proc(1) = wtest( y , h , kerntype , lagtype );
local t , w , sigma2 , e , x , b , m , kern , acov ;
    if (kerntype/=0) and (kerntype/=1);
       errorlog("wrong kernel choice");
       end;
    endif;
    T = rows(y);
    /* Calculating residuals from auxiliary regressions */
    if ( h == 0 );
       e = y;
    elseif ( h == 1 );
       x = ones(t,1);
       b = y/x;
       e = y-x*b;
    elseif ( h == 2 );
        x = ones(t,1)~seqa(1,1,t);
       b = y/x;
       e = y-x*b;
    else;
       errorlog("wrong null hypothesis: h");
    end;
    endif;
   /* Estimating autocovariances */
   acov = autocov(e);
   /* Calculating kernel */
   if ( lagtype[1] == 0 );
      m = lagtype[2];
   elseif ( lagtype[1] == 1 );
      m = lagsel( kerntype,acov,lagtype[2]);
   else;
      errorlog("lag selection variable: lagtype");
   end;
   endif;
   if (kerntype == 0);
       kern = bartlett( T , m );
   elseif ( kerntype == 1 );
       kern = qs( T , m );
   else;
      errorlog("wrong kernel selection: kerntype");
   end;
   endif;
   /* Estimating periodogram at frequency zero */
   m = 2.*ones(t,1);
   m[1] = 1;
   sigma2 = m'*(acov.*kern);
   /* Calculating test statistic */
   w = sumc(cumsumc(e)^2)./(T^2.*sigma2);
retp(w);
endp;
/*Procedure Bartlett:
** Purpose:Calculates Bartlett kernel.
** Format:kern = Bartlett( T , m )
** Input:
** T    scalar, sample size.
** m    scalar, bandwidth.
** Output:
** kern (Tx1)-vector, kernel weights for the estimated autocovariances for the estimation of
the smoothed periodogram.
** Note:
** k(m,j) = ( 1 - (j/(m+1)) )       for j <= m
**        = 0                       otherwise */
proc(1) = bartlett(T,m);
local j , kern;
     j = seqa(0,1,T);
     kern = ( 1 - j./(m+1));
     kern = kern.*(kern.>0);
retp(kern);
endp;

接上面。

/*** Procedure QS
** Purpose:Calculates Quadratic Spectral kernel.
** Format:kern = QS( T , m )
** Input:
** T    scalar, sample size.
** m    scalar, bandwidth.
** Output:
** kern (Tx1)-vector, kernel weights for the estimated autocovariances for the
**       estimation of the smoothed periodogram.
** Note:
** Quadratic spectral kernel in Andrews (1991,ECTA): "Heteroskedasticity and
** Autocorrelation Consistent Covariance Matrix Estimation". */
proc(1) = QS(T,m);
local j , kern;
    if ( m > 0 );
       j = seqa(1,1,T-1);
       j = j./m;
       kern = (25./(12.*pi^2.*j^2)).*( sin(1.2.*pi.*j)./(1.2.*pi.*j) - cos( 1.2.*pi.*j ) );
    elseif ( m == 0 );
       kern = zeros(T-1,1);
    else;
       errorlog("wrong bandwidth selection");
    end;
    endif;
    kern = 1|kern;
retp(kern);
endp;/*Procedure lagsel:
** Purpose:Automatic lag selection procedure from Newey and West (1994)
** (No prewhitening of the time series)
** Format:m = lagsel( kerntype , acov , n )
** Input:
** kerntype     scalar, = 0     Bartlett-kernel.
**                      = 1     Quadratic Spectral-kernel.
** acov         (Tx1)-vector, autocovariances of time series; output of
**                      procedure AUTOCOV.
** n            scalar, bandwidth parameter (integer).
** Output: m            scalar, bandwidth selected by procedure.
** Reference:
** Newey, W.K. and K.D. West (1994,RES), "Automatic Lag Selection in Covari-
** ance Matrix Estimation", pp. 631-653. */
proc(1) = lagsel( kerntype , acov , n );
local t , s0 , s1 , s2 , j , gam , m ;
    if ( kerntype /= 0 ) and ( kerntype /= 1 );
        errorlog("Wrong input of kernel choice");
    end;
    endif;
    if ( n > 0 );
       t = rows(acov);
       j = 2.*ones(n+1,1);
       j[1] = 1;
       s0 = j'*acov[1:n+1];
       /* Procedure for Bartlett kernel */
       if ( kerntype == 0 );
          j = 2.*seqa(0,1,n+1);
          s1 = j'*acov[1:n+1];
          gam = 1.1447.*((s1./s0)^2)^(1./3);
          m = minc(T|trunc(gam.*T^(1./3)));
       /* Procedure for Quadratic Spectral kernel */
       elseif ( kerntype == 1 );
          j = 2.*(seqa(0,1,n+1)^2);
          s2 = j'*acov[1:n+1];
          gam = 1.3221.*((s2./s0)^2)^(1./5);
          m = minc(T|(gam.*T^(1./5)));
       endif;
    elseif (n == 0 );
       m = 0;
    else;
       errorlog("wrong bandwidth choice");
    end;
    endif;
retp(m);
endp;

/*Procedure Nb:
** Purpose:
** Bandwidth sequence for Bartlett kernel with automatic lag selection.
** Sequence based on Newey and West (1994), table IIC, page 641.
** Format:n = nb(x,T);
** Input:
** x    scalar, parameter chosen by user.
** T    Sample size.
** Output:
** n    scalar, bandwith used as input for LAGSEL procedure. */
proc(1) = nb(x,T);
local n;
    n = trunc( x.*(T./100)^(2/9) );
retp(n);
endp;
/*Procedure Nqs:
** Purpose:
** Bandwidth sequence for Quadratic Spectral kernel with automatic
** lag selection.
** Sequence based on Newey and West (1994), table IIC, page 641.
** Format: n = nqs(x,T);
** Input:
** x    scalar, parameter chosen by user.
** T    Sample size.
** Output:
** n    scalar, bandwith used as input for LAGSEL procedure. */
proc(1) = nqs(x,T);
local n;
    n = trunc( x.*(T./100)^(2/25) );
retp(n);
endp;
/* Procedure Mb:
** Purpose:
** Bandwidth sequence for Bartlett kernel without automatic lag selection.
** Sequence based on Newey and West (1994), table IIC, page 641.
** Format: m = mb(x,T);
** Input:
** x    scalar, parameter chosen by user.
** T    Sample size.
** Output: m    scalar, bandwith used as input for WTEST procedure. */
proc(1) = mb(x,T);
local m;
      m = trunc( x.*(T./100)^(1/4) );
retp(m);
endp;
/* Procedure Mqs:
** Purpose:
** Bandwidth sequence for Quadratic spectral kernel without
** automatic lag selection. Sequence based on Hobijn, Franses and Ooms (1997).
** Format: m = mqs(x,T);
** Input:
** x    scalar, parameter chosen by user.
** T    Sample size.
Output: m    scalar, bandwith used as input for LAGSEL procedure. */
proc(1) = mqs(x,T);
local m;
    m = trunc( (2/3).*x.*(T./100)^(1/4) );
retp(m);
endp;
/* Procedure autocov:
** Purpose: Calculates autocovariances of time series.
** Format: gamk = autocov( y );
** Input: y    (Tx1)-vector, containing time series.
** Output: acov (Tx1)-vector, containing autocovariances of orders 0 through (T-1) */
proc(1) = autocov(y);
local k,t,yd,acov;
     t = rows(y);
     k = t-1;
     y = y - meanc(y);
     acov = rev(conv(y,rev(y),t-k,t));
retp(acov./T);
endp;

滞后项的选择非常重要，下面我们在考虑其他几种滞后项的选择方法，先从ADF检验的滞后项选择开始，然后在考虑KPSS检验，这些方法，一般都可以放在其他检验中滞后项的选择。

1、根据Schwert 提出的拇指准则选取， 设置lags=floor(12*(n/100)^(1/4))；
或者lags=floor(4*(n/100)^(1/4))；或者lags=floor(4*(n/100)^(2/9))。
还是用前面KPSS中的数据，下面的程序就省略了输入数据的语句，直接从调用函数开始，而且回归方程中带有漂移项，不带趋势项。

call adflags(y);
proc(0)=adfags(y);
local n,lags,dy0,y_1,dy,n1,k,x,tau,dy_i,i,label,R,q;
   n=rows(y);
   lags=floor(12*(n/100)^(1/4));    @取10个滞后@
   dy0=y[2:n]-y[1:n-1];
   y_1=y[1:n-1];
   /* lag 0 */
   dy=dy0[lags+1:n-1];
   y_1=y_1[lags+1:n-1];
   n1=rows(y_1); k=2;
   x=y_1~ones(n1,1);  @带趋势项  x=y_1~ones(n1,1)~seqa(1,1,n1); @
   tau=tOLS(dy,x);
   R={1 0,
      0 1};
   q={0,
      0};
   /* @如果带趋势项，估计方程是dy=delta*y_1+alpha+beta*t+Sum[gami*dy_i]+et.
    则约束条件用下面的语句。@
   R={1 0 0,
      0 1 0,
      0 0 1};
   q={0,
      0,
      0};
    */
   print/lz "n=" n1;
   label={"lags" "tau" "t dy_max" "p-value" "F stat"};
   print $label;
   print/rz 0~tau[1]~miss(0,0)~miss(0,0)~dfF(dy,x,R,q);
   /* lag 1 to lags */
   dy_i=ones(n1,1);
   i=1;                  
   do while i<=lags;                        @ 这个循环取10个滞后项的  @
      dy_i=dy_i~dy0[lags+1-i:n-1-i];
      x=y_1~dy_i;
      k=cols(x);
      tau=tOLS(dy,x);
      R=R~zeros(2,1);
      print/rz i~tau[1]~tau[k]~2*cdftc(abs(tau[k]),n1-k)~dfF(dy,x,R,q);   @2*cdftc(abs(tau[k]),n1-k) 是计算统计量的p值@
      i=i+1;
   endo;
  /* 我们把趋势项加入回归里。
   dy_i=ones(n1,1)~seqa(1,1,n1);         
   i=1;
   do while i<=lags;
      dy_i=dy_i~dy0[lags+1-i:n-1-i];    @dy0是y差分后的滞后@
      x=y_1~dy_i;
      k=cols(x);
      tau=tOLS(dy,x);
      print/rz i~tau[1]~tau[k]~2*cdftc(abs(tau[k]),n1-k);  
      i=i+1;
   endo;
*/


endp;

proc tOLS(y,x);     @ols回归调用函数，输出t值@
   local n1,k,b,e,s2hat,varb,se,t;
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'y;
   e=y-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   t=b./se;
   retp(t);
endp;
proc dfF(y,x,R,q);            @ols回归调用函数，输出F值@
   local n,k,nr,b,rssu,rssr,br,fstat;
   n=rows(x); k=cols(x); nr=rows(q);
   b=inv(x'x)*x'y;
   rssu=(y-x*b)'(y-x*b);
   br=b+inv(x'x)*R'inv(R*inv(x'x)*R')*(q-R*b);
   rssr=(y-x*br)'(y-x*br);
   fstat=((rssr-rssu)/nr)/(rssu/(n-k));
   retp(fstat);
endp;

结果为：
n= 39              
             lags              tau           t dy_max             p-value                 F stat
               0       -2.6544716                .                .                         4.2208164
               1       -2.7598758       0.91812825       0.36466239        4.4054904
               2       -2.4888546       -0.6546252       0.51698824        3.8651283
               3        -2.602916       0.91707459       0.36556321        4.0596264
               4       -2.5715032       0.42385468       0.67442197        4.0021604
               5       -2.4329103     0.0065598943       0.99480673        3.6876991
               6       -2.2166857      -0.26600835       0.79199404        3.2357239
               7       -2.2379659       0.47418382       0.63880209        3.2412406
               8       -1.9263461       -1.2376599       0.22576958        2.7036588
               9       -1.9301568      -0.98428148       0.33340647        2.9169332
              10       -1.9362932      -0.63686844        0.5295738        3.0291909

/*使用MAIC来确定最优的滞后长度
      sig2p=e'e/n1;
      taup=(b[1]^2)*(y_1'y_1)/sig2p;
      maic[i+1]=ln(sig2p)+2*(taup+i)/n1;
估计方程带趋势项和漂移项
*/
new; cls;
let data=
              80         10236000
              78         10351000
              71          9371000
              73          9126000
              70          8230000
              71         12111000
              72          5364000
              76          6222000
              76          3952000
              78          3859000
              80          4798000
              80          5133000
              77          4240000
              69         13687000
              69          5248000
              68          4058000
              67          5253000
              66          8963000
              67          6457000
              70          8882000
              69          4134000
              68          3736000
              69          3887000
              69          2095000
              68          1980000
              63          7433000
              63          5487000
              63          3996000
              58         10858000
              60         12095000
              57          7251000
              53         17632000
              55         42481000
              65         16051000
              64          9145000
              62         11264000
              64          4799000
              64          6888000
              66          7093000
              69          6598000
              70          6280000
              69          5372000
              65          7381000
              64          2792000
              63          4258000
              60          3074000
              62          4427000
              60          3767000
              61          3041000
              60         10240000
;
data=reshape(data,rows(data)/2,2);
y=data[.,1];
call adf3Ftestlags(y,10);

proc(0)=adf3Ftestlags(y,maxlags);
   local n,dy0,y_1,dy,n1,k,x,tau,dy_i,i,label,R,q,b,e,sig2p,taup,maic;
   maic=miss(zeros(maxlags+1,1),0);
   n=rows(y);
   dy0=y[2:n]-y[1:n-1];
   y_1=y[1:n-1];
   /* lag 0 */
   dy=dy0[maxlags+1:n-1];
   y_1=y_1[maxlags+1:n-1];
   n1=rows(y_1); k=2;
   x=y_1~ones(n1,1)~seqa(1,1,n1);
   tau=tOLS(dy,x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   sig2p=e'e/n1;
   taup=(b[1]^2)*(y_1'y_1)/sig2p;
   maic[1]=ln(sig2p)+2*(taup+0)/n1;
   R={1 0 0,
      0 1 0,
      0 0 1};
   q={0,
      0,
      0};
   print/lz "n=" n1;
   label={"lags" "tau" "t dy_max" "p-value" "F stat" "MAIC"};
   print $label;
   print/rz 0~tau[1]~miss(0,0)~miss(0,0)~dfF(dy,x,R,q)~maic[1];
   /* lag 1 to lags */
   dy_i=ones(n1,1)~seqa(1,1,n1);
   i=1;
   do while i<=maxlags;
      dy_i=dy_i~dy0[maxlags+1-i:n-1-i];
      x=y_1~dy_i;
      k=cols(x);
      tau=tOLS(dy,x);
      b=inv(x'x)*x'dy;
      e=dy-x*b;
      sig2p=e'e/n1; taup=(b[1]^2)*(y_1'y_1)/sig2p; maic[i+1]=ln(sig2p)+2*(taup+i)/n1;
      R=R~zeros(3,1);
      print/rz i~tau[1]~tau[k]~2*cdftc(abs(tau[k]),n1-k)~dfF(dy,x,R,q)~maic[i+1];
      i=i+1;
   endo;
   print/lz "# of opt.lags=" minindc(maic)-1;
endp;
proc tOLS(y,x);
   local n1,k,b,e,s2hat,varb,se,t;
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'y;
   e=y-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   t=b./se;
   retp(t);
endp;
proc dfF(y,x,R,q);
   local n,k,nr,b,rssu,rssr,br,fstat;
   n=rows(x); k=cols(x); nr=rows(q);
   b=inv(x'x)*x'y;
   rssu=(y-x*b)'(y-x*b);
   br=b+inv(x'x)*R'inv(R*inv(x'x)*R')*(q-R*b);
   rssr=(y-x*br)'(y-x*br);
   fstat=((rssr-rssu)/nr)/(rssu/(n-k));
   retp(fstat);
endp;
估计结果
n= 39              
              lags              tau            t dy_max             p-value                 F stat             MAIC
               0       -2.5410548                .                .                          2.9164461        72.656802
               1       -2.7577271         1.118239       0.27108169        3.2082349        94.450569
               2       -2.4078849      -0.43657408       0.66518027        2.7463572        84.388618
               3       -2.6566465        1.1637993        0.2528511        3.0870219        115.30698
               4       -2.6978027       0.69321158       0.49317788        3.1653202        143.99199
               5       -2.5875958       0.37427146       0.71075033        2.9961408        162.80864
               6       -2.3649908       0.09138603       0.92779335        2.6464045        168.05384
               7       -2.5336156       0.92381955       0.36320365        2.9000453        228.68146
               8        -1.906068       -0.7164363       0.47965544        2.1166835        170.13497
               9       -1.5355658      -0.61829336       0.54156313        2.0226486        132.12738
              10       -1.2660311      -0.40704143       0.68730923        1.9969107        107.42361
# of opt.lags= 0