# johansen test error in code

### Cross-posted at the Wolfram community forum

I have found this code on a forum, apparently it worked in a previous version of mathematica. Now when I run it, it keeps running for ages. I am not able to find the error.

 ewa = FinancialData["NYSE:EWA", {{2006, 4, 4}, {2012, 4, 9}}, "Value"];
ewc = FinancialData["NYSE:EWC", {{2006, 4, 4}, {2012, 4, 9}}, "Value"];
data = Transpose@{ewa, ewc};
johansenProcedure[levels_?MatrixQ, p_Integer] /;
(Dimensions[levels][[ 1]] > (Dimensions[levels][[2]] + p)) :=
With[{N = Dimensions[levels][[1]], T = Dimensions[levels][[2]],
diff = Rest[levels - RotateLeft[levels]], ypt = Drop[levels, -p]},
Module[{x, r0t, rpt, s00, sop, cc, λ, vλ, logλ},
x = Join @@ (Drop[RotateLeft[diff, #], p - 1] & /@ Range[p]);
With[{q = Inverse[Transpose[x].x], y0t = Drop[diff, p - 1]},
r0t = y0t - x.q.Transpose[x].y0t; rpt = ypt - x.q.Transpose[x].ypt; s00 = (Transpose[r0t].r0t)/T; s0p = (Transpose[r0t].rpt)/T;
cc = Inverse[CholeskyDecomposition[(Transpose[rpt].rpt)/T]];
{λ, vλ} = Eigensystem[(cc.Transpose[s0p].Inverse[s00].s0p.Transpose[cc])];
logλ = Log[1 - λ];
{Reverse[-T* FoldList[Plus, First[logλ], Rest[logλ]]],
Reverse[-T logλ], Join[Reverse[Transpose[Transpose[cc].vλ]], s00, s0p]} ]]];

johansenProcedure[data, 10]


This is the error code

Dot::dotsh: Tensors {<<1>>} and {{0.38,0.42},{-0.21,-0.27},{0.03,0.13},{0.,-0.03},<<44>>,      {0.04,0.02},{-0.56,-0.63},<<1455>>} have incompatible shapes. Dot::dotsh: Tensors {<<1>>} and {{0.38,0.42},{-0.21,-0.27},{0.03,0.13},{0.,-0.03},<<44>>,{0.04,0.02},{-0.56,-0.63},<<1455>>} have incompatible shapes.


On guy suggested that this line was broken:

r0t = y0t - x.q.Transpose[x].y0t;
rpt = ypt - x.q.Transpose[x].ypt;


I have spent ages but are not able to figure out what to do, could someone give it a try?

• There is a version of this code on @Verbeia 's website, using the old AppendRows function. I think this should be replaced with Join[##,2]& rather than just Join. Commented Sep 21, 2013 at 8:58
• Yes, and I already knew about the error from the Wolfram community forum. I haven't looked at this code for about 15 years. Due to work and family commitments, I'm not going to get around to fixing it for a while. Commented Sep 21, 2013 at 12:44
• Are there anyone els that are able to fix it? Commented Sep 21, 2013 at 13:54

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.

• Two questions 1) What is the source for the critical values 2) How does your code treat intercepts and trends? For example, there can be an intercept/ trend in the cointegrating vector or in the error correction model. I am asking so that I know which critical values should be used, e.g from Osterwald-Lenum (1992) (which are probably the same used in Stata). Yours are similar but not same to Panel 3. I can provide more detail if needed. Many thanks for posting the codes by the way, both to you and @Verbeia ! Commented Nov 10, 2018 at 11:30
• 1) I grabbed the critical values from the cointegration package source code jplv7 of the Econometrics Toolbox for MATLAB, spatial-econometrics.com. I haven't looked at those numbers for years. I should probably see if there is a better data source, such as the one you reference. Commented Nov 10, 2018 at 20:37
• 2) If I understand correctly, the optional deterministic term (intercept only) is in the error correction model. I append a constant offset to the xt array in the code, as you can see. (In the literature, this matrix is often labeled Z1.) However, I actually wrote this code a few years ago. I have an updated version that also allows for a trend term. This newer code can be found at: community.wolfram.com/groups/-/m/t/1397534 Commented Nov 10, 2018 at 20:48
• I should mention that at the time I added the critical values I wasn't aware that these values might depend upon the deterministic terms (intercept/slope) -- as apparently wasn't the author of the original paper. I should probably redo that part of the code with better values. However, it should be straightforward to modify it yourself, if you like. I should also mention that my Johansen code gives identical results to the Econometrics code. Commented Nov 10, 2018 at 20:58
• @Titus: I wrote code to generate the critical values, and I was able to reproduce the values (for the case of 3 time series, which is what I'm using in my portfolio) from my code (to within the error margin of the values). I did this for detrend = 0, 1, 2. (0 => No deterministic terms, 1 => unrestricted intercept, 2 => unrestricted intercept and trend.) So, I think they're probably good. Commented Nov 12, 2018 at 3:01

My interest piqued, I went back to my original code, from 2002 (ok, not quite 15 years). Simon's comment was correct. The problem in the modified version of the code (not mine) was that it changed the obsolete command AppendRows to Join instead of Join[##,2].

This version produces the expected format result in version 9, and does so in a fraction of a second. I also fixed up three other things - first, it no longer uses single capital letters for internal variables, which is bad practice; second, it replaces a FoldList with Accumulate, which was introduced in a version of Mathematica that post-dated the original version of this code. Thirdly and most importantly, I noticed another error which was that I was actually using the number of columns, not rows, of the matrix as t, the number of periods, so the result were wrong. The results below fix this. The version on my website is now updated: http://www.verbeia.com/mathematica/mma/johansenprocedure.nb

johansenProcedure[levels_?MatrixQ, p_Integer] /;
(Dimensions[levels][[1]] > (Dimensions[levels][[2]] + p)) :=
With[{n = Dimensions[levels][[2]], t = Dimensions[levels][[1]],
diff = Rest[levels - RotateLeft[levels]], ypt = Drop[levels, -p]},
Module[{x, r0t, rpt, s00, sop, cc, λ, vλ, logλ},
x = Join[##, 2] & @@ (Drop[RotateLeft[diff, #], p - 1] & /@  Range[p]);
With[{q = Inverse[Transpose[x].x], y0t = Drop[diff, p - 1]},
r0t = y0t - x.q.Transpose[x].y0t;
rpt = ypt - x.q.Transpose[x].ypt;
s00 = (Transpose[r0t].r0t)/t;
s0p = (Transpose[r0t].rpt)/t;
cc = Inverse[CholeskyDecomposition[(Transpose[rpt].rpt)/t]];
{λ, vλ} =
Eigensystem[(cc.Transpose[s0p].Inverse[s00].s0p.Transpose[cc])];
logλ = Log[1 - λ];
(* the result *)
{Reverse[-t*Accumulate[logλ]], Reverse[-t logλ],
Join[Reverse[Transpose[Transpose[cc].vλ]], s00, s0p]} ]]]


Testing:

  testdata = Table[Random[], {400}, {5}];
johansenProcedure[testdata, 2]
(*{{23.8744, 23.719, 23.5234, 21.1915, 13.039}, {0.155443, 0.195517,
2.33194, 8.15249,
13.039}, {{0.503254, -0.791089,
1.38498, -1.3733, -2.64845}, {-0.141704,
2.23088, -2.04149, -0.00648452, -1.87642}, {1.02837, -2.09406, \
-2.50085, 0.10001, -0.373228}, {-1.30376, 0.29492, 0.722075, 0.95909,
0.052498}, {0.13003, -0.777657, 0.356669,
3.08094, -1.28164}, {0.103774, -0.0086864, -0.00257725, \
-0.00449606, 0.0110278}, {-0.0086864, 0.107217, -0.00245552,
0.00366342, 0.00349774}, {-0.00257725, -0.00245552,
0.0935757, -0.000156411, -0.00914056}, {-0.00449606,
0.00366342, -0.000156411, 0.105339, -0.000924821}, {0.0110278,
0.00349774, -0.00914056, -0.000924821,
0.114438}, {-0.00450076, -0.00441737, 0.000598824,
0.0110404, -0.000333721}, {0.011656, 0.00527172,
0.0039008, -0.00258649, -0.000321094}, {0.00189628, -0.000695033,
0.00436633, -0.00120978, 0.00237625}, {0.00397597, -0.00905529,
0.00865361, -0.00282281, 0.00632063}, {-0.00523858,
0.000743599, -0.0087003, 0.000167249, 0.00202286}}} *)


Incidentally, the combination of a 2008-era computer instead of a 2000-era one, Mathematica version 9 not 4, fixing the bugs, and using Join instead of the obsolete AppendRows etc results in a speed-up of more than 100 times!

• Brilliant, thank you! You have no idea how much this helped me out! Commented Sep 22, 2013 at 10:56