This code is based on GAUSS code by Pedro de Lima of the Johns Hopkins University.
The input is the levels data, with each time period being one row (ie the matrix is t*n), and p, the number of lags in the VAR specification.
The output is a list containing:
 a vector containing the trace statistics to test H0: r<=q, H1: 0<=q<n
 a vector containing the lmax statistics to test H0: r<=q, H1: r=q+1
 a n*3 × n matrix containing the normalised eigenvectors, and the product moment matrices s00 and s0p, stacked one on top of each other.
 
 I have no idea how it's supposed to work, but it runs pretty fast (see Timing data below).

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-λ];
				{Reverse[-t*Accumulate[logλ]],Reverse[-t logλ],Join[Reverse[Transpose[Transpose[cc].vλ]],s00,s0p]} ]]]

testdata=Table[Random[],{400},{5}];Dimensions[testdata]

johansenProcedure[testdata,2]

This is the speed on my 2008-era computer, after updating to version-9 compatible code

Table[Timing[johansenProcedure[testdata,i];],{i,1,10}]

This was the speed on my 2002-era computer

testdata200=Table[Random[],{200},{5}];Dimensions[testdata200]

Table[Timing[johansenProcedure[testdata200,i];],{i,1,10}]

This was the speed on my 2002-era computer