new;

load xx[400,1] = Dwheatclose.out;
load zz[400,1] = Dcornclose.out;


xx = xx-meanc(xx);
zz = zz-meanc(zz);

y = xx~zz;

nobs=rows(y);


library optmum;

optset;


proc fct(a);
local ps, eps, phi, i, vc, corr, co, sco, epl, s1, pnobs, y11, y1l1, beps, ivc;

eps=zeros(nobs-1,4);

phi = zeros(2,2);

phi[1,1] = a[1];
phi[1,2] = a[2];
phi[2,1] = a[3];
phi[2,2] = a[4];


y11 = trimr(y,1,0); y1l1=trimr(y,0,1);
ps = y11' - phi * y1l1';


beps = ps'; 
eps[.,1:2] = beps;

eps[.,3] = (eps[.,1])^2;
eps[.,4] = (eps[.,2])^2;

eps[.,1] = eps[.,1]-meanc(eps[.,1]);
eps[.,2] = eps[.,2]-meanc(eps[.,2]);
eps[.,3] = eps[.,3]-meanc(eps[.,3]);
eps[.,4] = eps[.,4]-meanc(eps[.,4]);

corr = 0;

pnobs=rows(eps);

vc = vcx(eps);
ivc=inv(vc);


epl=eps;

i = 1;
do until i>10;

eps = trimr(eps,1,0);
epl=trimr(epl,0,1);


co = (eps'*epl)/(pnobs-1);

sco = co*ivc*co'*ivc;

corr = corr+ sumc(diag(sco));


i = i+1;
endo;

s1 = (pnobs)*(corr);





retp(s1);

endp;

start = { 0.1, 0.1, 0.5, 0.5 };


output file = cormaxw660ALT.out reset;

{ a, f, g, ret } = optprt(optmum(&fct, start));





end;