new;

load corn[300,1]= recorn.out;
load soy[300,1] = resoy.out;


nobs=200;
y = zeros(200,2);


i = 1;
do until i>(nobs-1);

t=nobs-i+1;
y[i,1] = corn[t];
y[i,2] = soy[t];
i = i+1;
endo;


output file= comod.out reset;
print y;

y = trimr(y, 0, 1);

library optmum;

optset;


proc fct(a);
local ps, eps, phi, i, term, vov, vc, corr, co, sco, epl, s1, s2, y11, y1l1, did, di, 
beps, nphi, y2l2, n1;

/*
eps= zeros(179,4);

*/

eps=zeros(197,6);

phi = zeros(2,2);
phi[1,1] = a[1];
phi[1,2] = a[2];
phi[2,1] = a[3];
phi[2,2] = a[4];

nphi = zeros(2,2);
nphi[1,1] = a[5];
nphi[1,2] = a[6];
nphi[2,1] = a[7];
nphi[2,2] = a[8];




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


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

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

eps[.,5]=(eps[.,3])^3; 
eps[.,6]=(eps[.,4])^3;


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]);
eps[.,5] = eps[.,5]-meanc(eps[.,5]);
eps[.,6] = eps[.,6]-meanc(eps[.,6]);




corr = zeros(6,6);

vc = vcx(eps);

did = diag(vc);
di = sqrt(did);
vov = di*di';


epl=eps;

i = 1;
do until i>10;

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

n1 = rows(eps);

co = (eps'*epl)/n1;

co = co./vov;

sco = co^2;

corr = corr+sco;


i = i+1;
endo;

s1 = sumc(corr);


s2 = sumc(s1);





retp(s2);

endp;


/*
start = { 0.5, 0.5, -0.1, 0.1, 0.5, 0.5, -0.1, 0.1 };
*/


start = { 1.3, -0.2, 1.3, 0.6, -0.1, 0.1, -1.2, 0.1 };

output file = VARmax2.out reset;

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


hh=_opfhess;
v=inv(hh);
se=sqrt(diag(v));
tt=a./se;

q=a~se~tt;


print q;

end;