options crt;

?  autoregressive order is 2
?  Generates Johansen and Juselius results with Finnish data.
?  originally by David Cushman
?  VAR version   by Clint Cummins 11/1/93
?  corrections to eigenvectors, and normalization  1/28/94

load;   ? load data from bottom of program

dummy;
? To create "centered" seasonal dummies
? Not actually necessary unless no "linear trend in the nonstationary part of
? the data generating process," as in the Johansen and Juselius model for the
? Danish data.

d1=q1-q4; d2=q2-q4; d3=q3-q4;

set t=@nob;

print t;

set n=4;

dot 1-4; 
smpl 58:3 84:3;
genr dy.=y.-y.(-1);
smpl 58:4 84:3;
genr ly.=y.(-2);
enddot;

? The VAR command is an easy way to regress different dependent variables
? on the same independent variables and generate residual covariance matrices.

VAR(silent) dy1-dy4 | c dy1(-1) dy2(-1) dy3(-1) dy4(-1) d1 d2 d3;
rename @res veo;
mat S00 = @covu*[(@nob-8)/@nob];   ? divide veo'veo by @NOB instead of @NOB-K

VAR(silent) ly1-ly4 | c dy1(-1) dy2(-1) dy3(-1) dy4(-1) d1 d2 d3;
rename @res vek;
mat Skk = @covu*[(@nob-8)/@nob];

set t1=@nob;
print t1;

mat S0k=(veo'vek)/t1;

?  Johansen-Juselius Theorem 4.1, eqn (4.5) is:
?  |lam*Skk - Sk0*S00"*S0k| = 0
?  This determinantal equation can be converted into a standard eigenvalue
?  problem for a symmetric matrix, which TSP can handle, by pre- and post-
?  multiplying by the inverse matrix square root of Skk, denoted ccinv below:
?  cc = chol(Skk), Skk = cc'cc,  ccinv = cc",  Skk" = ccinv*ccinv'
?  so  ccinv'Skk*ccinv = ccinv'*cc'cc*ccinv = I*I = I  , and
?  |lam*I - ccinv'*Sk0*S00"*S0k*ccinv| = 0
?  so  lam  is the eigenvalues of  ccinv'*Sk0*S00"*S0k*ccinv .
?  Note that the eigenvectors need to be transformed by ccinv.

?  The eigenvalue problem in Theorem 5.1 can be handled with the same type
?  of pre- and post- multiplication by the inverse square root of H'Skk*H .

?  The matrix algebra below combines several steps from the original JOH1.TSP
?  program.  JOH1 used several steps due to a bug in TSP's MAT command in
?  early 4.2A versions (those dated before 5/23/91).

mat ccinv = (chol(skk))";
mat B=ccinv'(S0k'S00"S0k)*ccinv;
print B;
?  Get the eigenvalues of B (real and positive) in descending order.
mat vec=ccinv*eigvec(B);   ? transform for eigenvectors of original problem
rename @eigval eigen; ? this is a way in 4.2B to avoid doing the eigenvalue/
                      ? eigenvector calculation twice.  In 4.2A or earlier, use:
? mat eigen=eigval(B);
print eigen vec;  ? see also below for sorted and normalized versions

freq n;
smpl 1 n;
unmake eigen eigs;
sort(reverse) eigs;
print eigs;   ? print sorted eigenvalues

?  Normalize the eigenvectors so they look like J-J Table 2

mat nd = ident(n);
set nd(1,1) = 1/vec(1,1);
set nd(2,2) = 1/vec(4,2);
set nd(3,3) = 1/vec(1,3);
set nd(4,4) = 1/vec(1,4);
mat vecnormd = vec*nd;
print vecnormd;

?  Construct likelihood ratio test statistics
ll=log(1-eigs);
set j=1;
dot 0-3;
 smpl j n;
 msd(noprint) ll; 
 set lr. = -t1*@sum;
 set j=j+1;
enddot;

print lr0 lr1 lr2 lr3;

end;

?  Finnish data for Johansen and Juselius results.

noprint;
freq q;
smpl 58:2 84:3;
load y1;     ? real money
         2.848664 2.828010 2.908186
2.934326 2.986127 3.006543 3.016235
2.954030 3.002750 3.041873 3.033884
3.058968 3.035220 3.087247 3.100822
3.058580 3.039497 3.075370 3.107436
3.137721 3.127759 3.128404 3.199932
3.068053 3.072799 3.084895 3.135006
3.099091 3.090548 3.078614 3.110883
3.024568 3.083357 3.041080 3.117137
3.026565 3.046901 2.993750 3.006729
2.991557 3.047116 3.073535 3.180678
3.135102 3.193558 3.196119 3.333479
3.411379 3.459963 3.451505 3.530157
3.297392 3.280233 3.307619 3.451498
3.437284 3.485342 3.500827 3.606705
3.508857 3.580815 3.494112 3.668480
3.495950 3.553317 3.536008 3.689291
3.646950 3.663690 3.565684 3.783132
3.641322 3.636200 3.580232 3.647325
3.629197 3.645420 3.585797 3.595339
3.599353 3.693067 3.663110 3.693763
3.661466 3.789403 3.752624 3.813215
3.725020 3.770843 3.713289 3.745108
3.722089 3.733176 3.732510 3.788271
3.705715 3.813328 3.821703 3.849595
3.803960 3.863014 3.852309 3.826078
3.782735 3.830368 3.825409;

load y2;       ? real income
         3.836539 3.902982 3.896880
3.847045 3.916454 3.986216 3.984211
3.956458 3.990079 4.048979 4.053306
4.066436 4.088541 4.146216 4.125302
4.108145 4.095878 4.157463 4.159411
4.100841 4.169761 4.199699 4.197109
4.178529 4.178418 4.192901 4.260383
4.222265 4.219907 4.276832 4.293151
4.199678 4.276213 4.299289 4.341603
4.293609 4.315433 4.322150 4.312939
4.305840 4.333206 4.380844 4.392300
4.420402 4.444100 4.489679 4.518776
4.490783 4.536882 4.561972 4.619116
4.510130 4.575117 4.569593 4.631627
4.602710 4.623280 4.643740 4.715684
4.689491 4.691851 4.703741 4.809920
4.745103 4.809414 4.751507 4.880777
4.762205 4.816440 4.757928 4.807621
4.695385 4.772790 4.766532 4.842663
4.695384 4.732926 4.729782 4.801802
4.698796 4.762504 4.739455 4.843240
4.771751 4.841930 4.816520 4.890468
4.785139 4.837455 4.867240 4.922232
4.799556 4.875043 4.868998 4.936803
4.802846 4.890625 4.891206 4.968430
4.834531 4.923077 4.935910 4.980688
4.866376 4.959666 4.952320;

load y3;     ? marginal interest rate
        .173113 .101202 .132080
.069992 .065319 .065319 .065319
.065319 .161268 .243495 .163309
.110109 .122483 .117338 .147644
.175297 .229126 .242397 .240355
.210909 .171261 .189132 .170671
.136714 .168645 .177979 .197210
.183321 .135405 .067659 .067659
.067659 .115380 .220019 .294459
.067659 .121864 .091850 .120889
.091667 .067659 .067659 .067659
.067659 .067659 .067659 .071948
.082317 .151347 .231112 .185400
.079273 .117961 .092853 .081580
.074644 .074644 .074644 .074644
.074644 .074644 .095583 .231112
.168476 .147816 .089932 .137760
.129887 .175465 .161949 .246235
.178565 .177979 .159394 .167462
.141326 .152464 .165514 .180570
.170755 .091485 .087095 .092853
.080566 .080381 .081119 .111899
.116093 .139066 .148334 .153150
.163309 .121775 .124428 .140370
.130941 .129360 .128217 .141065
.133919 .139762 .144533 .155122
.161268 .154436 .151691;

load y4;      ? inflation rate
        .014815 .000000 .000000
.000000 .000000 .014599 .014389
.014185 .000000 .013986 .000000
.000000 .000000 .000000 .027399
.000000 .026668 .000000 .013072
.012903 .012739 .012579 .012423
.059898 .022990 .011300 .000000
.011173 .021979 .000000 .010811
.010695 .010582 .020834 .010257
.010152 .010050 .009950 .029270
.037740 .018349 .000000 .009050
.000000 .008969 .000000 .000000
.017700 .000000 .008734 .008658
.025533 .024898 .024293 .007968
.007905 .031010 .015152 .014926
.029199 .035339 .054067 .025975
.043894 .036149 .057487 .016621
.058684 .035627 .044017 .028304
.045462 .013245 .038715 .020878
.040491 .035091 .026467 .011132
.014652 .018019 .014185 .007018
.027588 .020203 .016529 .019481
.034759 .042560 .026433 .025752
.030600 .032348 .018397 .012937
.032873 .024571 .007255 .021455
.019386 .037235 .012623 .027351
.004274 .019008 .014538;