wfopen C:\Users\amanor\Downloads\stokk.wf1

' set sample 
sample s0 2/01/2008 2/10/2018
sample s1 2/05/2008 2/10/2018


' initialization of parameters and starting values
' change below only to change the specification of model 
smpl s0


'get starting values from univariate GARCH 
equation eq1.ARCH(CGARCH, OPTMETHOD=LEGACY, MAXIT=1000, fixedlag=32, MAXIT=1000) Y1 C AR(1)
equation eq2.ARCH(1, 2, CGARCH, OPTMETHOD=LEGACY, fixedlag=32) Y2 C AR(1)

' declare coef vectors to use in bi-variate GARCH model
' see above for details 

'C-GARCH(1,1): GARCH = Q + alpha(RESID(-1)^2 - Q(-1)) + beta(GARCH(-1) - Q(-1)), Q = omega + rho(Q(-1) - omega) + phi(RESID(-1)^2 - GARCH(-1)) can be rewritten as a GARCH(2,2) process:

'GARCH(2,2) = omega + alpha(1)*RESID(-1)^2 + alpha(2)*RESID(-2)^2 + beta(1)*GARCH(-1) + beta(2)*GARCH(-2) 

'Hence: 
'C-GARCH(1,1) = GARCH(2,2) = (1-alpha-beta)(1-rho)*omega + (alpha + phi)*RESID(-1)^2 - (alpha*rho + (alpha + beta)*phi)*RESID(-2)^2 + (beta - phi)*GARCH(-1) - (beta*rho - (alpha + beta)*phi)*GARCH(-2)

'omega(1) = (1-alpha-beta)(1-rho)*omega
'alpha(1) = (alpha + phi) - (alpha*rho +(alpha + beta)*phi)
'beta(1) = (beta - phi) - (beta*rho - (alpha + beta)*phi)

coef(2) mu
mu(1) = eq1.c(1)
mu(2)= eq2.c(1)

coef(3) omega
omega(1)= ((1-eq1.c(6)-eq1.c(7))*(1-eq1.c(4))*eq1.c(3))^.5
omega(2)=0.01
omega(3)=((1-eq2.c(6)-eq2.c(7))*(1-eq2.c(4))*eq2.c(3))^.5

coef(4) alpha
alpha(1) = @abs(eq1.c(6) + eq1.c(5) - (eq1.c(6)*eq1.c(4) + ((eq1.c(6) + eq1.c(7))*eq1.c(5))))^.5
alpha(2) = 0
alpha(3) = 0
alpha(4) = @abs(eq2.c(6) + eq2.c(5) - (eq2.c(6)*eq2.c(4) + ((eq2.c(6) + eq2.c(7))*eq2.c(5))))^.5

coef(4) beta 
beta(1)= @abs((eq1.c(7) - eq1.c(5)) - (eq1.c(7)*eq1.c(4) - (eq1.c(6)+eq1.c(7))*eq1.c(5)))^.5
beta(2)= 0
beta(3) = 0
beta(4)= @abs((eq2.c(7) - eq2.c(5)) - (eq2.c(7)*eq2.c(4) - (eq2.c(6)+eq2.c(7))*eq2.c(5)))^.5

' constant adjustment for log likelihood
!mlog2pi = 2*log(2*@acos(-1))

' use var-cov of sample in "s1" as starting value of variance-covariance matrix
series cov_y1y2 = @cov(y1-mu(1), y2-mu(2))
series var_y1 = @var(y1)
series var_y2 = @var(y2)

series sqres1 = (y1-mu(1))^2
series sqres2 = (y2-mu(2))^2
series res1res2 = (y1-mu(1))*(y2-mu(2))


' ...........................................................
' LOG LIKELIHOOD
' set up the likelihood 
' 1) open a new blank likelihood object (L.O.) name bvgarch
' 2) specify the log likelihood model by append
' ...........................................................

logl bvgarch
bvgarch.append @logl logl
bvgarch.append sqres1 = (y1-mu(1))^2
bvgarch.append sqres2 = (y2-mu(2))^2
bvgarch.append res1res2 = (y1-mu(1))*(y2-mu(2))

' calculate the variance and covariance series

bvgarch.append var_y1 = omega(1)^2 + alpha(1)^2*sqres1(-1)+ 2*(alpha(1)*alpha(3)*res1res2(-1))+alpha(3)^2*sqres2(-1)+beta(1)^2*var_y1(-1)+2*(beta(1)*beta(3)*cov_y1y2(-1))+beta(3)^2*var_y2(-1)

bvgarch.append var_y2 = omega(2)^2+omega(3)^2 + alpha(2)^2*sqres2(-1)+2*(alpha(2)*alpha(4)*res1res2(-1))+alpha(4)^2*sqres1(-1)+beta(2)^2*var_y1(-1)+2*(beta(2)*beta(4)*cov_y1y2(-1))+beta(4)^2*var_y2(-1)

bvgarch.append cov_y1y2 = omega(2)*omega(1) + alpha(1)*alpha(2)*sqres1(-1) + (alpha(2)*alpha(3) + alpha(1)*alpha(4))*res1res2(-1) + alpha(3)*alpha(4)*sqres2(-1) +beta(1)*beta(2)*var_y1(-1) + (beta(2)*beta(3) + beta(1)*beta(4))*cov_y1y2(-1) + beta(3)*beta(4)*var_y2(-1) 

' determinant of the variance-covariance matrix
bvgarch.append deth = var_y1*var_y2 - cov_y1y2^2

' inverse elements of the variance-covariance matrix
bvgarch.append invh1 = var_y2/deth
bvgarch.append invh3 = var_y1/deth
bvgarch.append invh2 = -cov_y1y2/deth

' log-likelihood series
bvgarch.append logl =-0.5*(!mlog2pi + (invh1*sqres1+2*invh2*res1res2+invh3*sqres2) + log(deth))

' remove some of the intermediary series
bvgarch.append @temp invh1 invh2 invh3 sqres1 sqres2 res1res2 deth

' estimate the model
smpl s1
bvgarch.ml(showopts, optmethod=legacy, m=2000, c=0.01)

' change below to display different output
show bvgarch.output
graph varcov.line var_y1 var_y2 cov_y1y2
show varcov

' LR statistic for univariate versus bivariate model
scalar lr = -2*( eq1.@logl + eq2.@logl-bvgarch.@logl)
scalar lr_pval = 1 - @cchisq(lr,1)