I've also written a Mathematica version of the Johansen test which I coded after studying the equations for the procedure. I've examined Johansen test codes (in Matlab, C, etc.) available online. I've found that these different codes often don't agree with each other, especially with respect to detrending and normalization.
My code allows for detrending of the data (0 => no detrend, 1 => remove offset). I usually use offset detrending. In addition to detrending, I've included a short routine to display the Johansen statistics in a nice readable form. The input is m x N data (which is the transpose of the verbeia.com code), where m is the number of vectors, and N is the length of the vectors. The output is an array containing an array of eigenvalues, an array of eigenvectors, and the statistics array.
JohansenTest[data_List, order_Integer, detrend_Integer] :=
Module[{i, k, diff, cc, uRes, vRes, tLen, suu, svv, suv, svu, val, vec, jMat, num, len, lenDiff, tmp, lr1, lr2, cvm, cvt, stats, xt = {}, yt, debug = False, sgns, vecs, ord},
If[order < 2,
Print["Johansen Test error: order must be >= 2 in order to have lagged difference terms. Aborting..."]; Abort[]];
If[(detrend < 0) || (detrend > 1),
Print["JohansenRegress Error: Invalid detrend parameter, ", detrend, ". Must be 0 <= detrend <= 1. Aborting..."; Abort[]]];
(* For the VAR(order) model, k = order-1 is number of lagged difference terms in corresponding VECM *)
k = order - 1;
{num, len} = Dimensions[data];
diff = data[[All, k + 1 ;;]] - data[[All, k ;; -2]];
lenDiff = Dimensions[diff][[2]];
Do[xt = Join[xt, diff[[All, k - i + 1 ;; lenDiff - i]]], {i, 1, k}];
If[detrend > 0, AppendTo[xt, ConstantArray[1, {lenDiff - k}]]];
xt = xt\[Transpose];
(* Regress dy(t) against dy(t-1), dy(t-2), ..., dy(t-k) *)
yt = diff[[All, k + 1 ;;]]\[Transpose];
uRes = yt - xt.Inverse[xt\[Transpose].xt].xt\[Transpose].yt;
(* Regress y(t-k) against dy(t-1), dy(t-2),..., dy(t-k) *)
yt = data[[All, 2 ;; -k - 1]]\[Transpose];
vRes = yt - xt.Inverse[xt\[Transpose].xt].xt\[Transpose].yt;
tLen = Length[uRes];
svv = vRes\[Transpose].vRes/tLen;
suu = uRes\[Transpose].uRes/tLen;
suv = uRes\[Transpose].vRes/tLen;
svu = suv\[Transpose];
(* Calculate Cholesky decomposition of svv. Used to pre- and post-multiply eiqenvector equation *)
(* LowerTriangularize function zeroes out any (small) non-zero elements in the upper triangle (due rounding errors from the Inverse calc). *)
cc = LowerTriangularize[Inverse[CholeskyDecomposition[svv]\[Transpose]]];
If[debug, Print["Normalization check 1: ", Chop[cc.svv.cc\[Transpose]]//MatrixForm]];
jMat = cc.svu.Inverse[suu].suv.cc\[Transpose];
(* The eigenvectors from 'Eigensystem' are rows, not columns: *)
{val, vec} = Eigensystem[jMat];
(* Normalize eigenvectors so that vec.svv.vec\[Transpose] = identity matrix *)
vec=(cc\[Transpose].vec)\[Transpose];
If[debug, Print["Normalization check 2: ",
Chop[vec.svv.vec\[Transpose]] // MatrixForm]];
If[debug, Print["Unit Root Test p-value for ", num, " eigenvector-weighted data: ", {Range[num], UnitRootTest[vec[[#]].data] & /@ Range[num]}\[Transpose] // TableForm]];
(* Sort eigenvalues and vectors from largest to smallest eigenvalue *)
vec = vec[[Reverse[Ordering[val]]]];
val = Reverse[Sort[val]];
(* Compute statistics *)
lr1 = Table[0, {num}];
lr2 = Table[0, {num}];
cvm = Table[0, {num}, {3}];
cvt = Table[0, {num}, {3}];
Do[
tmp = Log[1 - val][[i ;;]];
lr1[[i]] = -tLen*(Plus @@ tmp);
lr2[[i]] = -tLen*Log[1 - val[[i]]];
cvm[[i]] = csja[num - i + 1, -1];
cvt[[i]] = csjt[num - i + 1, -1];
, {i, 1, num}
];
stats = {lr1, lr2, cvt, cvm};
{val, vec, stats}
]
PrintStats[stats_List]:=
Module[{null},
null=Table["r<="<>ToString[i-1],{i,1,Length[stats[[1]]]}];
Print["NULL\tTrace Statistic\tCrit 90%\tCrit 95%\tCrit 99%\n",null//TableForm,"\t",stats[[3,1]]//TableForm,"\t\t\t",stats[[3,3]]//TableForm];
Print["NULL\tEigen Statistic\tCrit 90%\tCrit 95%\tCrit 99%\n",null//TableForm,"\t",stats[[3,2]]//TableForm,"\t\t\t",stats[[3,4]]//TableForm];
]
csja[n_Integer,p_Integer]:=
Module[{jcp,out},
jcp[0]={{2.9762`,4.1296`,6.9406`},{9.4748`,11.2246`,15.0923`},{15.7175`,17.7961`,22.2519`},{21.837`,24.1592`,29.0609`},{27.916`,30.4428`,35.7359`},{33.9271`,36.6301`,42.2333`},{39.9085`,42.7679`,48.6606`},{45.893`,48.8795`,55.0335`},{51.8528`,54.9629`,61.3449`},{57.7954`,61.0404`,67.6415`},{63.7248`,67.0756`,73.8856`},{69.6513`,73.0946`,80.0937`}};
jcp[1]={{2.7055`,3.8415`,6.6349`},{12.2971`,14.2639`,18.52`},{18.8928`,21.1314`,25.865`},{25.1236`,27.5858`,32.7172`},{31.2379`,33.8777`,39.3693`},{37.2786`,40.0763`,45.8662`},{43.2947`,46.2299`,52.3069`},{49.2855`,52.3622`,58.6634`},{55.2412`,58.4332`,64.996`},{61.2041`,64.504`,71.2525`},{67.1307`,70.5392`,77.4877`},{73.0563`,76.5734`,83.7105`}};
jcp[2]={{2.7055`,3.8415`,6.6349`},{15.0006`,17.1481`,21.7465`},{21.8731`,24.2522`,29.2631`},{28.2398`,30.8151`,36.193`},{34.4202`,37.1646`,42.8612`},{40.5244`,43.4183`,49.4095`},{46.5583`,49.5875`,55.8171`},{52.5858`,55.7302`,62.1741`},{58.5316`,61.8051`,68.503`},{64.5292`,67.904`,74.7434`},{70.463`,73.9355`,81.0678`},{76.4081`,79.9878`,87.2395`}};
If[(p>1)||(p<-1)||(n>12)||(n<1),
out=Table[0,{3}],
out=jcp[p+1][[n]]
];
out
]
csjt[n_Integer,p_Integer]:=
Module[{jcp,out},
jcp[0]={{2.9762`,4.1296`,6.9406`},{10.4741`,12.3212`,16.364`},{21.7781`,24.2761`,29.5147`},{37.0339`,40.1749`,46.5716`},{56.2839`,60.0627`,67.6367`},{79.5329`,83.9383`,92.7136`},{106.7351`,111.7797`,121.7375`},{137.9954`,143.6691`,154.7977`},{173.2292`,179.5199`,191.8122`},{212.4721`,219.4051`,232.8291`},{255.6732`,263.2603`,277.9962`},{302.9054`,311.1288`,326.9716`}};
jcp[1]={{2.7055`,3.8415`,6.6349`},{13.4294`,15.4943`,19.9349`},{27.0669`,29.7961`,35.4628`},{44.4929`,47.8545`,54.6815`},{65.8202`,69.8189`,77.8202`},{91.109`,95.7542`,104.9637`},{120.3673`,125.6185`,135.9825`},{153.6341`,159.529`,171.0905`},{190.8714`,197.3772`,210.0366`},{232.103`,239.2468`,253.2526`},{277.374`,285.1402`,300.2821`},{326.5354`,334.9795`,351.215`}};
jcp[2]={{2.7055`,3.8415`,6.6349`},{16.1619`,18.3985`,23.1485`},{32.0645`,35.0116`,41.0815`},{51.6492`,55.2459`,62.5202`},{75.1027`,79.3422`,87.7748`},{102.4674`,107.3429`,116.9829`},{133.7852`,139.278`,150.0778`},{169.0618`,175.1584`,187.1891`},{208.3582`,215.1268`,228.2226`},{251.6293`,259.0267`,273.3838`},{298.8836`,306.8988`,322.4264`},{350.1125`,358.719`,375.3203`}};
If[(p>1)||(p<-1)||(n>12)||(n<1),
out=Table[0,{3}],
out=jcp[p+1][[n]]
];
out
]
The main code is JohansenTest. It calls the PrintStats routine displays the statistics, for example:
I welcome feedback. This code gives different results than the verbeia.com code, even when I use no detrending, and I don't think the difference is only a matter of normalization. If I have time, I'll compare the two and see what is going on.
EDIT: I think I see the main difference in the two codes. The way I read the equations, the line in the verbeia.com code
x = Join[##, 2] & @@ (Drop[RotateLeft[diff, #], p - 1] & /@ Range[p]);
should be
x = Join[##, 2] & @@ (Drop[RotateLeft[diff, #], p - 1] & /@ Range[p-1]);
where Range[p] has been replaced with Range[p-1]. If the author is still around, perhaps he could comment on this. Am I correct? This is what the equations I have show. This would also bring his code into agreement with my code (up to that point) except for an overall factor of -1 for the matrices and the fact that my matrices are the transpose of his -- both of which probably don't matter.
EDIT #2: I forgot to include the two small functions csja and csjt, so I just added them. Apologies.
EDIT #3: As of 8/15/2017, I've rewritten the code to simplify it and also correct an eigenvector normalization issue. I've also added a few debug lines that can be helpful to make sure that the code is working correctly. Specifically, the "normalization check" lines should both display an identity matrix, which insures that that the eigenvector is correctly normalized. (The correctly normalized eigenvectors are the rows of the 'vec' matrix.) Finally, the last debug line displays the p-value of a unit root test applied to the eigenvector-weighted data, for each of the eigenvectors. These p-values will indicate the cointegration rank. For example, p-values <= 0.01 indicate a stationary (cointegrated) array with 99% confidence. If, say, two of the eigenvectors produce low p-values, but the other eigenvector(s) do not, then this indicates two cointegrating relationships.
AppendRows
function. I think this should be replaced withJoin[##,2]&
rather than justJoin
. $\endgroup$