load corn-oil.wf1

'set series
series y1= 100*dlog(oil)
series y2= 100*dlog(corn)


'declare coefficients
coef(2) beta0
beta0(1)=0.0
beta0(2)=0.0

coef(4) beta
beta(1)=0.0
beta(2)=0.0
beta(3)=0.0
beta(4)=0.0

coef(2) lambda
lambda(1)=0.2
lambda(2)=0.2

coef(2) alpha0
alpha0(1)=0.001
alpha0(2)=0.001

coef(2) gamma
gamma(1)=0.5
gamma(2)=0.5

coef(4) alpha
alpha(1)=0.001
alpha(2)=0.3
alpha(3)=0.001
alpha(4)=0.2

coef(2) delta
delta(1)=-0.01
delta(2)=-0.01

coef(2) ef
ef(1)=0.2
ef(2)=0.1

coef(1) c0
c0(1)=0.01

coef(1) ci
ci(1)=0.04

coef(1) g
g(1)=0.7

!pi = @acos(-1)
!absz=sqr(2/!pi)

'set presample values of expressions in logl

series res1=y1-beta0(1)-beta(1)*y1(-1)-beta(2)*y2
series res2=y2-beta0(2)-(beta(3))*y1(-1)-beta(4)*y2(-1)
!var1=(@stdev(res1))^2
!var2=(@stdev(res2))^2
!correlation=@cor(y1,y2)

series h11=!var1
series h22=!var2
series rho=!correlation
series h12=2*rho*sqr(h11)*sqr(h22)
series deth=h11*h22-h12*h12
series invh11=h22/deth
series invh12=-h12/deth
series invh22=h11/deth


'set up EGARCH likelihood

logl ll1
ll1.append @logl logl
@byeqn
ll1.append res1=y1-beta0(1)-beta(1)*y1(-1)-beta(2)*y2
ll1.append res2=y2-beta0(2)-beta(3)*y1(-1)-beta(4)*y2(-1)

ll1.append h11=exp(alpha0(1)+gamma(1)*log(h11(-1))+alpha(1)*(abs(res1(-1)/sqr(h11(-1)))-!absz+delta(1)*((res1(-1))/sqr(h11(-1)))+alpha(2)*(abs(res2/sqr(h22))-!absz+delta(2)*(res2/sqr(h22)))

ll1.append h22=exp(alpha0(2)+gamma(2)*log(h22(-1))+alpha(3)*(abs(res2(-1)/sqr(h22(-1)))-!absz+delta(2)*((res2(-1)/sqr(h22(-1)))+(alpha(4))*(abs(res1(-1)/sqr(h11(-1)))-!absz+delta(1)*((res1(-1))/sqr(h11(-1))))

ll1.append rho=@cor(res1,res2)
ll1.append h12=rho*sqr(h11)*sqr(h22)
ll1.append deth=h11*h22-h12*h12
ll1.append invh11=h22/deth
ll1.append invh12=-h12/deth
ll1.append invh22=h11/deth

ll1.append likelihoods=log(deth)+(res1^2)*invh11+(res2^2)*invh22+2*res1*res2*invh12
ll1.append logl=-log(2*!pi)-0.5*likelihoods

ll1.append z1=res1/sqr(h11)
ll1.append z2=res2/sqr(h22)
ll1.append z11=(res1^2)/h11
ll1.append z22=(res2^2)/h22
ll1.append z12=z1*z2

' estimate and display output

ll1.ml(showopts,m=100,c=1e-5)
show ll1.output