' bv_garch_m 

' model
'mean equation 1        :    y1 = c1(1)*(h1,t)^0.5 + c2(1)*(h2,t)^0.5 + c3(1) + res1
'mean equation 2        :    y2 = c1(2)*(h1,t)^0.5 + c2(2)*(h2,t)^0.5 + c3(2) + res2
'variance equation 1   :    h1,t = c4(1) + c5(1)*(res1,t-1)^2 + c6(1)*h1,t-1
'variance equation 2   :    h2,t = c4(2) + c5(2)*(res2,t-1)^2 + c6(2)*h2,t-1

'THE AIM IS TO CAPTURE THE EFFECTS OF (h1,t)^0.5 AND (h2,t)^0.5 IN THE MEAN EQUATION.

'change path to program path
%path = @runpath
cd %path

'load workfile
'wfload intl_fin.wf1

wfopen "bv_data.wf1"

'dependent variables of both series must be continues
smpl @all
series y1 = dp_mex 'The name of the first series
series y2 = dq_mex 'The name of the second series

wfsave "bv_garch_m_results.wf1"

' sample first observation of s1 need to be one or two periods after the first observation of s0 
sample s0 @first @last-5 'MONTH-1 MONTH-T
sample s1 @first+1 @last-5 'MONTH-2 MONTH-T

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

'get starting values from univariate GARCH-M
equation eq1.arch(archm=sd, conv=1e-5, maxit=100, fixedlag=32, conv=1e-5, maxit=100) y1 c
equation eq2.arch(archm=sd, conv=1e-5, maxit=100, fixedlag=32, conv=1e-5, maxit=100) y2 c

'save the conditional variances
eq1.makegarch h1
eq2.makegarch h2

'declare coef vectors to use in bi-variate GARCH model
coef(2) c1 '@sqrt(h1) term
c1(1) = eq1.c(1) 
c1(2) = 0

coef(2) c2 '@sqrt(h2) term
c2(1) = 0
c2(2) = eq2.c(1)

coef(2) c3 'c term in the mean equation
c3(1) = eq1.c(1) 
c3(2) = eq2.c(1)

coef(2) c4 ' c term in the garch equation
c4(1) = eq1.c(3)
c4(2) = eq2.c(3)

coef(2) c5 'resid(-1)^2 term in the garch equation
c5(1) = eq1.c(4)
c5(2) = eq2.c(4)

coef(2) c6 'garch(-1) term in the garch equation
c6(1) = eq1.c(5)
c6(2) = eq2.c(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 var_y1 = @var(y1)
series var_y2 = @var(y2)
series cov_y1y2 = @cov(y1-c1(1)*h1^0.5-c3(1), y2-c2(2)*h2^0.5-c3(2))

series res1 = y1-c1(1)*h1^0.5-c3(1)
series res2 = y2-c2(2)*h2^0.5-c3(2)
series res1res2 = (y1-c1(1)*h1^0.5-c3(1))*(y2-c2(2)*h2^0.5-c3(2))

'SPECIFY LOG LIKELIHOOD
logl bv_garch_m
bv_garch_m.append @logl logl '@byeqn '@byobs
bv_garch_m.append res1 = y1-c1(1)*h1^0.5-c2(1)*h2^0.5-c3(1)
bv_garch_m.append res2 = y2-c1(2)*h1^0.5-c2(2)*h2^0.5-c3(2)
bv_garch_m.append res1res2 = (y1-c1(1)*h1^0.5-c2(1)*h2^0.5-c3(1))*(y2-c1(2)*h1^0.5-c2(2)*h2^0.5-c3(2))
bv_garch_m.append var_y1 = c4(1) + c5(1)*res1^2 + c6(1)*var_y1(-1)
bv_garch_m.append var_y2 = c4(2) + c5(2)*res2^2 + c6(2)*var_y2(-1)
bv_garch_m.append cov_y1y2 = (c4(1) + c5(1)*res1^2 + c6(1)*var_y1(-1)) * (c4(2) + c5(2)*res2^2 + c6(2)*var_y2(-1))

'determinant of the variance-covariance matrix and inverse elements of the variance-covariance matrix
bv_garch_m.append deth = var_y1*var_y2 - cov_y1y2^2
bv_garch_m.append invh1 = var_y1/deth
bv_garch_m.append invh2 = var_y2/deth
bv_garch_m.append invh12 = -cov_y1y2/deth

' log-likelihood series
bv_garch_m.append logl = -0.5*(!mlog2pi + (invh2*(res1^2)+2*invh12*(res1res2)+invh1*(res2^2)) + log(deth)) '"buraya bak"

'estimate the model
'smpl s2
'var_y1=1
'var_y2=1

smpl s1
bv_garch_m.ml(showopts, m=100, c=1e-5) '"maximum likelihood tahmin et"

group chk var_y1 var_y2 cov_y1y2 res1 res2 deth invh1 invh12 invh2

'change below to display different output
show bv_garch_m.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 - bv_garch_m.@logl )
scalar lr_pval = 1 - @cchisq(lr,1)