'change path to program path
 %path = @runpath
 cd %path
 load macro.wf1

 ' dependent variables of both series must be continues
 smpl @all
 series y1 = DSPXM
 series y2 = DFTSE100M
 series y3 = DIPUSM
 series y4 = DIPUKM

 ' set sample
 ' first observation of s1 need to be one or two periods after
 ' the first observation of s0
sample s2 1990M01 2014M12
sample s3 1990M02 2014M12

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

 'get starting values from univariate GARCH
 equation eq1.arch(m=100,c=1e-5) y1 c
 equation eq2.arch(m=100,c=1e-5) y2 c
 equation eq3.arch(m=100,c=1e-5) y3 c
 equation eq4.arch(m=100,c=1e-5) y4 c

 ' declare coef vectors to use in bi-variate GARCH model
 ' see above for details
 coef(4) mu
 mu(1)= eq1.c(1)
 mu(2)= eq2.c(1)
 mu(3)= eq3.c(1)
 mu(4)= eq4.c(1)

 coef(4) teta
 teta(1) = .1
 teta(2)= .1
 teta(3)= .1
 teta(4)= .1

 coef(10) omega
 omega(1)=eq1.c(3)^0.5
 omega(2)=0.1
 omega(3)=0.1
 omega(4)=0.1

 omega(5)=eq2.c(3)^0.5
 omega(6)=0.1
 omega(7)=0.1

 omega(8)=eq3.c(3)^0.5
 omega(9)=0.1

 omega(10)=eq3.c(3)^0.5

 coef(4) alpha
 alpha(1) = 0.1
 alpha(2) = 0.2
 alpha(3) = 0.3
 alpha(4) = 0.3
 
 coef(4) beta
 beta(1) = 0.3
 beta(2) = 0.3
 beta(3) = 0.3
 beta(4) = 0.3

 ' constant adjustment for log likelihood
 !mlog2pi = 4*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 cov_y1y3 = @cov(y1-mu(1), y3-mu(3))
 series cov_y1y4 = @cov(y1-mu(1), y4-mu(4))
 series cov_y2y3 = @cov(y2-mu(2), y3-mu(3))
 series cov_y2y4 = @cov(y2-mu(2), y4-mu(4))
 series cov_y3y4 = @cov(y3-mu(3), y4-mu(4))
 series var_y1 = @var(y1)
 series var_y2 = @var(y2)
 series var_y3 = @var(y3)
 series var_y4 = @var(y4)

 series sqres1 = (y1-mu(1))^2
 series sqres2 = (y2-mu(2))^2
 series sqres3 = (y3-mu(3))^2
 series sqres4 = (y4-mu(4))^2

 series res1res2 = (y1-mu(1))*(y2-mu(2))
 series res1res3 = (y1-mu(1))*(y3-mu(3))
 series res1res4 = (y1-mu(1))*(y4-mu(4))
 series res2res3 = (y2-mu(2))*(y3-mu(3))
 series res2res4 = (y2-mu(2))*(y4-mu(4))
 series res3res4 = (y3-mu(3))*(y4-mu(4))

 smpl 1990M01 1990M02

 series res1=y1-mu(1)

 series res2=y2-mu(2)

 series res3=y3-mu(3)

 series res4=y4-mu(4)

 smpl 1990M02 2014M04

 series res1=y1-mu(1)-teta(1)*res1(-1)
 series res2=y2-mu(2)-teta(2)*res2(-1)
 series res3=y3-mu(3)-teta(3)*res3(-1)
 series res4=y4-mu(4)-teta(4)*res4(-1)

 series cov_y1y2 = @cov(y1-mu(1)-teta(1)*res1(-1), y2-mu(2)-teta(2)*res2(-1))
 series cov_y1y3 = @cov(y1-mu(1)-teta(1)*res1(-1), y3-mu(3)-teta(3)*res3(-1))
 series cov_y1y4 = @cov(y1-mu(1)-teta(1)*res1(-1), y4-mu(4)-teta(4)*res4(-1)) 
 series cov_y2y3 = @cov(y2-mu(2)-teta(2)*res2(-1), y3-mu(3)-teta(3)*res3(-1))
 series cov_y2y4 = @cov(y2-mu(2)-teta(2)*res2(-1), y4-mu(4)-teta(4)*res4(-1))
 series cov_y3y4 = @cov(y3-mu(3)-teta(3)*res3(-1), y4-mu(4)-teta(4)*res4(-1))

 series var_y1 = @var(y1-teta(1)*res1(-1))
 series var_y2 = @var(y2-teta(2)*res2(-1))
 series var_y3 = @var(y3-teta(3)*res3(-1))
 series var_y4 = @var(y4-teta(4)*res4(-1))
 series sqres1=y1-mu(1)-teta(1)*res1(-1)^2
 series sqres2=y2-mu(2)-teta(2)*res2(-1)^2
 series sqres3=y3-mu(3)-teta(3)*res3(-1)^2
 series sqres4=y4-mu(4)-teta(4)*res4(-1)^2
 series res1res2=(y1-mu(1)-teta(1)*res1(-1))*(y2-mu(2)-teta(2)*res2(-1))
 series res1res3=(y1-mu(1)-teta(1)*res1(-1))*(y3-mu(3)-teta(3)*res3(-1))
 series res1res4=(y1-mu(1)-teta(1)*res1(-1))*(y4-mu(4)-teta(4)*res4(-1))
 series res2res3=(y2-mu(2)-teta(2)*res2(-1))*(y3-mu(3)-teta(3)*res3(-1))
 series res2res4=(y2-mu(2)-teta(2)*res2(-1))*(y4-mu(4)-teta(4)*res4(-1))
 series res3res4=(y3-mu(3)-teta(3)*res3(-1))*(y4-mu(4)-teta(4)*res4(-1))

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

 logl fvgarch
 fvgarch.append @logl logl

 smpl 1990M01 1990M02
 fvgarch.append res1=y1-mu(1)
 fvgarch.append res2=y2-mu(2) 
 fvgarch.append res3=y3-mu(3) 
 fvgarch.append res4=y4-mu(4) 

 smpl 1990M01 1990M02

 fvgarch.append sqres1=(y1-mu(1))^2
 fvgarch.append sqres2=(y2-mu(2))^2
 fvgarch.append sqres3=(y3-mu(3))^2
 fvgarch.append sqres4=(y4-mu(4))^2
 fvgarch.append res1res2= (y1-mu(1))*(y2-mu(2))
 fvgarch.append res1res3= (y1-mu(1))*(y3-mu(3))
 fvgarch.append res1res4= (y1-mu(1))*(y4-mu(4))
 fvgarch.append res2res3= (y2-mu(2))*(y3-mu(3))
 fvgarch.append res2res4= (y2-mu(2))*(y4-mu(4))
 fvgarch.append res3res4= (y3-mu(3))*(y4-mu(4))


 smpl 1990M02 2014M12

 fvgarch.append res1=y1-mu(1)-teta(1)*res1(-1)
 fvgarch.append res2=y2-mu(2)-teta(2)*res2(-1)
 fvgarch.append res3=y3-mu(3)-teta(3)*res3(-1)
 fvgarch.append res4=y4-mu(4)-teta(4)*res4(-1)


 smpl 1990M02 2014M12

 fvgarch.append sqres1 = (y1-mu(1)-teta(1)*res1(-1))^2
 fvgarch.append sqres2 = (y2-mu(2)-teta(2)*res2(-1))^2
 fvgarch.append sqres3 = (y3-mu(3)-teta(3)*res3(-1))^2
 fvgarch.append sqres4 = (y4-mu(4)-teta(4)*res4(-1))^2
 fvgarch.append res1res2 = (y1-mu(1)-teta(1)*res1(-1))*(y2-mu(2)-teta(2)*res2(-1))
 fvgarch.append res1res3 = (y1-mu(1)-teta(1)*res1(-1))*(y3-mu(3)-teta(3)*res3(-1))
 fvgarch.append res1res4 = (y1-mu(1)-teta(1)*res1(-1))*(y4-mu(4)-teta(4)*res4(-1))
 fvgarch.append res2res3 = (y2-mu(2)-teta(2)*res2(-1))*(y3-mu(3)-teta(3)*res3(-1))
 fvgarch.append res2res4 = (y2-mu(2)-teta(2)*res2(-1))*(y4-mu(4)-teta(4)*res4(-1))
 fvgarch.append res3res4 = (y3-mu(3)-teta(3)*res3(-1))*(y4-mu(4)-teta(4)*res4(-1))

 ' calculate the variance and covariance series
fvgarch.append var_y1 = omega(1)^2 + beta(1)^2*var_y1(-1) + alpha(1)^2*sqres1(-1)
fvgarch.append var_y2 = omega(2)^2+omega(5)^2 + beta(2)^2*var_y2(-1) + alpha(2)^2*sqres2(-1)
fvgarch.append var_y3 = omega(3)^2+omega(6)^2+omega(8)^2+ beta(3)^2*var_y3(-1) + alpha(3)^2*sqres3(-1)
fvgarch.append var_y4 = omega(4)^2+omega(7)^2+omega(9)^2+omega(10)^2+ beta(4)^2*var_y4(-1) + alpha(4)^2*sqres4(-1)

fvgarch.append cov_y1y2 = omega(1)*omega(2)* + beta(2)*beta(1)*cov_y1y2(-1) + alpha(2)*alpha(1)*res1res2(-1)
fvgarch.append cov_y1y3 = omega(1)*omega(3) + beta(3)*beta(1)*cov_y1y3(-1) + alpha(3)*alpha(1)*res1res3(-1)
fvgarch.append cov_y1y4 = omega(1)*omega(4) + beta(4)*beta(1)*cov_y1y4(-1) + alpha(4)*alpha(1)*res1res4(-1)
fvgarch.append cov_y2y3 = omega(2)*omega(3) + omega(4)*omega(5) + beta(3)*beta(2)*cov_y2y3(-1) + alpha(3)*alpha(2)*res2res3(-1)
fvgarch.append cov_y2y4 = omega(2)*omega(4) + omega(5)*omega(6) + beta(4)*beta(2)*cov_y2y4(-1) + alpha(4)*alpha(2)*res2res4(-1)
fvgarch.append cov_y3y4 = omega(3)*omega(5) + omega(6)*omega(7) + beta(3)*beta(4)*cov_y3y4(-1) + alpha(3)*alpha(4)*res3res4(-1)

 ' determinant of the variance-covariance matrix
 fvgarch.append deth = (var_y1*((var_y2*((var_y3*var_y4)-(cov_y3y4^2)))-(cov_y2y3*((cov_y2y3*var_y4)-(cov_y3y4*cov_y2y4)))+(cov_y2y4*((cov_y2y3*cov_y3y4)-(var_y3*cov_y2y4)))))-(cov_y1y2*((cov_y1y2*((var_y3*var_y4)-(cov_y3y4^2)))-(cov_y2y3*((cov_y1y3*var_y4)-(cov_y3y4*cov_y1y4)))+(cov_y2y4*((cov_y1y3*cov_y3y4)-(var_y3*cov_y1y4)))))+(cov_y1y3*((cov_y1y2*((cov_y2y3*var_y4)-(cov_y3y4*cov_y2y4)))-(var_y2*((cov_y1y3*var_y4)-(cov_y3y4*cov_y1y4)))+(cov_y2y4*((cov_y1y3*cov_y2y4)-(cov_y2y3*cov_y1y4)))))-(cov_y1y4*((cov_y1y2((cov_y2y3*cov_y3y4)-(var_y3*cov_y2y4)))-(var_y2*((cov_y1y3*cov_y3y4)-( var_y3*cov_y1y4)))+(cov_y2y3*((cov_y1y3*cov_y2y4)-(cov_y2y3*cov_y1y4)))))

 ' inverse elements of the variance-covariance matrix
 fvgarch.append invh1 = (cov_y2y3*cov_y3y4*cov_y2y4 - cov_y2y4*var_y3*cov_y2y4 + cov_y2y4*cov_y2y3*cov_y3y4 - var_y2*cov_y3y4*cov_y3y4 - cov_y2y3*cov_y2y3*var_y4 + var_y2*var_y3*var_y4)/deth
 fvgarch.append invh2 = (cov_y2y4*var_y3*cov_y1y4 - cov_y2y3*cov_y3y4*cov_y1y4 - cov_y2y4*cov_y1y3*cov_y3y4 + cov_y1y2*cov_y3y4*cov_y3y4 + cov_y2y3*cov_y1y3*var_y4 - cov_y1y2*var_y3*var_y4)/deth
 fvgarch.append invh3 = (var_y2*cov_y3y4*cov_y1y4 - cov_y2y4*cov_y2y3*cov_y1y4 + cov_y2y4*cov_y1y3*cov_y2y4 - cov_y1y2*cov_y3y4*cov_y2y4 - var_y2*cov_y1y3*var_y4 + cov_y1y2*cov_y2y3*var_y4)/deth
 fvgarch.append invh4 = (cov_y2y3*cov_y2y3*cov_y1y4 - var_y2*var_y3*cov_y1y4 - cov_y2y3*cov_y1y3*cov_y2y4 + cov_y1y2*var_y3*cov_y2y4 + var_y2*cov_y1y3*cov_y3y4 - cov_y1y2*cov_y2y3*cov_y3y4)/deth
 fvgarch.append invh5 = (cov_y1y3*cov_y3y4*cov_y1y4 - cov_y1y4*var_y3*cov_y1y4 + cov_y1y4*cov_y1y3*cov_y3y4 - var_y1*cov_y3y4*cov_y3y4 - cov_y1y3*cov_y1y3*var_y4 + var_y1*var_y3*var_y4)/deth
 fvgarch.append invh6 = (cov_y1y4*cov_y2y3*cov_y1y4 - cov_y1y2*cov_y3y4*cov_y1y4 - cov_y1y4*cov_y1y3*cov_y2y4 + var_y1*cov_y3y4*cov_y2y4 + cov_y1y2*cov_y1y3*var_y4 - var_y1*cov_y2y3*var_y4)/deth
 fvgarch.append invh7 = (cov_y1y2*var_y3*cov_y1y4 - cov_y1y3*cov_y2y3*cov_y1y4 + cov_y1y3*cov_y1y3*cov_y2y4 - var_y1*var_y3*cov_y2y4 - cov_y1y2*cov_y1y3*cov_y3y4 + var_y1*cov_y2y3*cov_y3y4)/deth
 fvgarch.append invh8 = (cov_y1y2*cov_y2y4*cov_y1y4 - cov_y1y4*var_y2*cov_y1y4 + cov_y1y4*cov_y1y2*cov_y2y4 - var_y1*cov_y2y4*cov_y2y4 - cov_y1y2*cov_y1y2*var_y4 + var_y1*var_y2*var_y4)/deth
 fvgarch.append invh9 = (cov_y1y3*var_y2*cov_y1y4 - cov_y1y2*cov_y2y3*cov_y1y4 - cov_y1y3*cov_y1y2*cov_y2y4 + var_y1*cov_y2y3*cov_y2y4 + cov_y1y2*cov_y1y2*cov_y3y4 - var_y1*var_y2*cov_y3y4)/deth
 fvgarch.append invh10 = (cov_y1y2*cov_y2y3*cov_y1y3 - cov_y1y3*var_y2*cov_y1y3 + cov_y1y3*cov_y1y2*cov_y2y3 - var_y1*cov_y2y3*cov_y2y3 - cov_y1y2*cov_y1y2*var_y3 + var_y1*var_y2*var_y3)/deth

 ' log-likelihood series
 fvgarch.append logl =-0.5*(!mlog2pi + (invh1*sqres1+2*invh2*res1res2+2*invh3*res1res3+2*invh4*res1res4+invh5*sqres2+2*invh6*res2res3+2*invh7*res2res4+ invh8*sqres3+2*invh9*res3res4+invh10*sqres4) + log(deth)) '

' remove some of the intermediary series
 ' fvgarch.append @temp invh1 invh2 invh3 sqres1 sqres2 res1res2 deth
 ' estimate the model
 smpl s3
 fvgarch.ml(showopts, m=500, c=1e-5)

 ' change below to display different output
 show fvgarch.output
 graph var.line var_y1 var_y2 var_y3 var_y4
 graph cov.line cov_y1y2 cov_y1y3 cov_y1y4 cov_y2y3 cov_y2y4 cov_y3y4
 show var
 show cov

 ' LR statistic for univariate versus fourvariate model
 scalar lr = -2*(eq1.@logl + eq2.@logl + eq3.@logl + eq4.@logl - fvgarch.@logl )
 scalar lr_pval = 1 - @cchisq(lr,4)