new;

& simulates a bivariate noncausal VAR(1) with t-distributed errors @
& and estimates with gcov

@ parameters @
/*
a = { 0.3 0.02, 0.1 1.1 };
b = inv(a);
*/

a = { 1 -1, 0 1 };
b=inv(a);

a1 = a[.,1];
a2 = a[.,2];






b1 = b[1,.];
b2 = b[2,.];





jey1 = 0.7;

jey2 = 2;

jej = jey1|jey2;
ej = eye(2).*jej;
/*
nobs=400;
*/

nobs=200;

@ Simulation starts here @


x= zeros(nobs,2);

nu1 = 4;
nu2= nu1;

z1 = rndn(nobs,2);

i = 1;
do until i >nobs;

z2 = rndn(nu1,1);
sz2 = z2^2;
suz= sumc(sz2);
nsu = suz/nu1;
sq1 = sqrt(nsu);
if sq1 == 0;
goto here;
endif;


az2 = rndn(nu2,1);
asz2 = az2^2;
asuz= sumc(asz2);
ansu = asuz/nu2;
sq2 = sqrt(ansu);

if sq2 == 0;
goto here;
endif;

x[i,1] = z1[i,1]/sq1;
x[i,2] = z1[i,2]/sq2;


here:
i = i+1;

endo;


y1 = zeros(nobs,1);

y1[1] = b1*x[1,.]';

i = 2;
do until i>nobs;
y1[i] = jey1*y1[i-1] + b1*x[i,.]';
i = i+1;
endo;

y2 = zeros(nobs,1);

y2[nobs] = b2*x[nobs,.]';

i = 1;

do until i> (nobs-1);
y2[nobs-i] = (1/jey2)*y2[nobs-i+1] - (1/jey2)*b2*x[nobs-i+1,.]';

i = i+1;

endo;


y = zeros(nobs,2);

i = 1;
do until i>nobs;
opa = a1*y1[i,1] + a2*y2[i,1];
y[i,.] = opa';
i = i+1;
endo;

@ Estimation starts here @

library optmum;

optset;

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


te=nobs-1;
eps=zeros(te,4);
phi = zeros(2,2);
phi[1,1] = d[1];
phi[1,2] = d[2];
phi[2,1] = d[3];
phi[2,2] = d[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 = zeros(4,4);

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);


co = (eps'*epl)/nobs;

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 };
*/

start = { 0.7, -1.3, 0.0, 2.0 };

output file = cormax.out reset;

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


output file = sim.out reset;

print d';

pdi=zeros(2,2);

pdi[1,1]=  d[1];
pdi[1,2]=  d[2];
pdi[2,1]=  d[3];
pdi[2,2]=  d[4];


{ va, ve } = eigv(pdi);



print va';



end;