/***  Program for Monte Carlo examination of Johansen procedure ***/

   library gauss coint;
   clear dx, xk, z, q, r0t, rkt, sig,a,b,d,e,k,lr1,lr2, zz ;
   clear s00, sk0, skk, tmp, crit,i,f ;
   clear one, z_s, xk_s, c, cinv, a1, a2, d1, d2, lr3  ;
   clear m00,m01,m0k,m10,m11,m1k,mk0,mk1,mkk,t ;

/**/


/**/
output file = johan.out reset ;
/* Generate series */
/**/

capT = 400 ;
np=2 ;                 @dimension of VAR@
 y = zeros(capT,np) ;
 temp_y = zeros(capT,np) ;
seed1=91697 ;
 temp_y[1:3,.] = rndns(3,np,seed1) ; @initial values for process - here random@


let sigmamat= {.1 0 , 0 .1} ;  @variance matrix for innovations@
let pimat=  {0 0 , 0 0} ;  @The Pi matrix in Johansen notation@
PI1 = Pimat + eye(np) ;   @I+Pimat to simulat y_t = (I+Pi)*y_{t-1} +...@ 
sigsqrt = chol(sigmamat) ;
sigsqrt = sigsqrt' ;       @gauss creates upper triangular. Prefer lowr tr@
mu = zeros(2,1) ; @0.1*ones(2,1) ;@

@critical values for Johansen test. See Oxford Bull. August '92, p. 468@
crval95 = {277.71, 233.13, 192.89, 156.00, 124.24, 94.15,
           68.53, 47.21, 29.68, 15.41, 3.76} ;
crval99 = {293.44, 247.18, 204.95, 168.36, 133.57, 103.18,
           76.07, 54.46, 35.65, 20.04, 6.65} ;

p=-1 ;  @parameter for adjusting the data. Used in the SJstar procedure@
        @-1:unchanged;   0:demean;   1:detrend; .... @

level = 95 ;  @testing at 1% or 5%@
Nmc = 50 ;    @Number of estimations in Monto Carlo @
imc = 0 ;
Nuroots = zeros(Nmc,3) ;
Nu_k    = zeros(Nmc,1) ;
Nu_ka   = zeros(Nmc,1) ;
bic_v = zeros(3,2);
aic_v = zeros(3,2);

/****/
if level == 95 ;
      crval = crval95 ;
      "level = " ;; level ;
elseif level == 99 ;
      crval = crval99 ;
      "level = " ;; level ;
endif ;
/****/

"true number of unit roots in estimated model: " ;; 2 ;
" In DGP the no. of lags used is k=1 " ;
datestr(date) ;
timestr(time) ;

do while imc < Nmc ;    @beginning of Monte Carlo loop @

  imc = imc +1 ;
  p_count = imc % 100 ;
  if p_count == 0.0 ;
     imc ;             @to print out how many replications every 100th one @
  endif ;
  ii = 3 ;
  do while ii < capT ;
   ii = ii + 1 ;
seed=91697+capt*imc+ii-1 ;
   eps = sigsqrt*rndns(np,1,seed) ;                 @error term@
   temp_y[ii,.] = mu' + temp_y[ii-1,.]*PI1' + eps';
  endo ;
  y[.,.] = temp_y[1:capT,.] ;
  nn=cols(y) ;

  kk=1 ;
  do while kk<4 ;
   k=kk ;
   format /rd 2,0 ;
   format /mb1 /ROS 16,8 ;
   {eval, evec, evecs, lr1, lr2,  cvt, cvm} = sjstar(y,p,k,1) ;


/*    /* Calc. AIC and BIC from SJ residuals */
   nparam = nn + k*nn*nn ;
   temp = ln(abs(det(S00))) + sumc(ln(1-eval)) ;
   aic = temp + 2*nparam/rows(y) ;
   bic = temp + nparam*ln(rows(y))/rows(y) ;
      bic_v[kk,2] = bic ;
   aic_v[kk,2] = aic ;               */

   jj=1 ;              @ lr1[1] is the largest eigenvalue @
   no_rt = nn ;
   do until jj > nn ;
     if ( lr1[jj] < crval[11-nn+jj] ) ;
       goto done ;
     else ;
       no_rt = no_rt-1 ;
       jj=jj+1 ;
     endif ;
   endo ;
   done:
   nuroots[imc,kk] = no_rt ;

   kk=kk+1 ;
  endo ;  /* end of kk loop */


endo ; /* end of monte carlo loop ****/


@After Nmc replications add up the # of times 1,2,..or, nn roots are found@
xxx= ones(nmc,1) ;

kkk=1;
do while kkk < 4;
  " " ;
  " k = " kkk ;
  p_r=0 ;
  do until p_r > nn ;
   "frequency of" p_r "uroots: " ;; sumc(nuroots[.,kkk] .== p_r*xxx) ;
   p_r = p_r + 1 ;
  endo ;
  kkk = kkk+1 ;
endo ;

" " ;
timestr(time) ;
output off ;


/* SJstar.ARC
   to do Johansen  Econometrica 11/1991 test with constant mu
   (see also Johansen and Juselius 1990 Oxford Bulletin)
   H1: PI has rank p and mu.ne.0
   H2: PI=alpha beta'   has rank r<p
   H2*: PI=alpha beta'   has rank r<p   and    mu=alpha beta * */

proc (7) = sjstar(x,p,k,r);

/* x is a Tx1 matrix.  k=# of lags.
   z denotes the matrix containing the constant term
   z_s denotes the z* matrix;  as in the paper  */

/*     if (p > -1) ;    @to remove a trend - not used @
        f = 0 ;
     else ;
        f = p ;
     endif ;
     x     = detrend(x,p);     */
     dx    = diff(x,1) ;
     if k > 1 ;
        q     = 2 ;                 @k must be at least 2 to make sense@
        z     = lag(dx,1) ;
        do until (q >= k);
           z  = z~lag(dx,q);
           q  = q + 1 ;
        endo ;
        z_s = z ;
        one = ones(rows(z),1) ;
        z = z_s~one ;
     else ;
        z_s = zeros(rows(dx),1) ;
        one = ones(rows(dx),1) ;
        z= one ;
     endif ;
     f = -1 ;  @ no demeaning or detrending here @

     z     = detrend(trimr(z,k,0),f);
     z_s     = detrend(trimr(z_s,k,0),f);
     dx    = detrend(trimr(dx,k,0),f);

     xk    = detrend(trimr(lag(x,k),k+1,0),f);
     one = ones(rows(xk),1) ;
     xk_s = xk~one ;

     t = rows(xk_s) ;
     m00 = dx'dx/t ;
     m01 = dx'z/t ;
     m0k = dx'xk/t ;
     m10 = z'dx/t ;
     m11 = z'z/t ;
     m1k = z'xk/t ;
     mk0 = xk'dx/t ;
     mk1 = xk'z/t ;
     mkk = xk'xk/t ;
     r0t   = dx - (m01*inv(m11)*z')' ; @ resids of DX on DX(lagged),1 @
     rkt   = xk - (mk1*inv(m11)*z')' ; @ resids of X(-k) on DX(lagged),1 @
     skk   = rkt'rkt/rows(rkt) ;
     sk0   = rkt'r0t/rows(rkt) ;
     s00   = r0t'r0t/rows(r0t) ;
     sig   = sk0*inv(s00)*(sk0');
     c = chol(skk) ; c = c' ;  @transpose since GAUSS gives the transpose@
     cinv = inv(c) ;
     {a,d} = eigrs2(cinv*sig*cinv'/1);

     /* Normalize the eigen vectors such that (d'skk*d) = I */
     d     = cinv'd ;

/* Order the eigen values and eigen vectors together */

     zz    = (-a)~(d') ;          @negative, to sort in decreasing order @
     zz    = sortc(zz,1) ;
     a1     = -zz[.,1] ;
     d1     = zz[.,2:cols(zz)]' ;  /* Done */

/* Compute the trace statistics
      -2Q = -T SUM{1 to r} log[1-lambda(i)] r=1,..p  */

     lr1   =  zeros(rows(a1),1);    @LR1[ r=p , r=p-1 , ....r=2, r=1 ] @
     i     =  1 ;
   @"a1 " ;; a1 ; pause(2) ;@
   do while (i <= rows(a1));
        lr1[i,1] = -rows(rkt)*sumc(trimr(ln(1-a1),i-1,0));
        i = i + 1 ;
     endo ;

/*   For the hypo of H(r|r+1), the test stat is:   */
     lr2   = -rows(rkt)*ln(1-a1[r+1:r+1,1]) ;


retp(-sortc(-a,1),d1,d2,lr1,lr2,c_sjt(cols(x)-r,p),c_sja(cols(x)-r,p));
endp ;


end;