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. –  Simon Woods Sep 21 '13 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. –  Verbeia Sep 21 '13 at 12:44
Are there anyone els that are able to fix it? –  ALEXANDER Sep 21 '13 at 13:54

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! –  ALEXANDER Sep 22 '13 at 10:56