' y = mu + res -> y = H*lambda + mu + res
' res ~ N(0,H)
' H = omega*omega' + beta H(-1) beta' + alpha res(-1) res(-1)' alpha'
' where y = 2 x 1
' mu = 2 x 1
' lambda = 2 x 1
' H = 2 x 2 (symmetric)
' H(1,1) = variance of y1 (saved as var_y1)
' H(1,2) = cov of y1 and y2 (saved as cov_y1y2)
' H(2,2) = variance of y2 (saved as var_y2)
' omega = 2 x 2 low triangular
' beta = 2 x 2 diagonal
' alpha = 2 x 2 diagonal

smpl @all
series y1 = p1
series y2 = pm
' set sample for GARCH estimation (not the whole series > leave some observations for forecasting)
sample s0 1 2247 
sample s1 2 2247
' load data
smpl s0
'get starting values for parameters from univariate GARCH-M (1,1); archm=var shows the inclusion of
var in the mean eq
equation eq1.arch(archm=var,m=100,c=1e-5) y1 c
equation eq2.arch(archm=var,m=100,c=1e-5) y2 c
'save the conditional variances
eq1.makegarch garch1
eq2.makegarch garch2
' declare coef vectors to use in bi-variate GARCH model (please see introduction for details)
coef(2) lambda
lambda(1) = eq1.c(1)
lambda(2) = eq2.c(1)
coef(2) mu
mu(1) = eq1.c(2)
mu(2)= eq2.c(2)
coef(3) omega
omega(1)=(eq1.c(3))^.5
omega(2)=0 ' because we don't have it in the univariate GARCH-M
omega(3)=(eq2.c(3))^.5
coef(2) alpha
alpha(1) = (eq1.c(4))^.5
alpha(2) = (eq2.c(4))^.5
coef(2) beta
beta(1)= (eq1.c(5))^.5
beta(2)= (eq2.c(5))^.5
' constant adjustment for log likelihood (i.e. we define 2log(2pi))
!mlog2pi = 2*log(2*@acos(-1))
'old values
' 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))
series cov_y1y2 = @cov(y1-mu(1)-lambda(1)*garch1, y2-mu(2)-lambda(2)*garch2)
series var_y1 = @var(y1-lambda(1)*garch1)
series var_y2 = @var(y2-lambda(2)*garch2)
series sqres1 = (y1-mu(1)-lambda(1)*garch1)^2
series sqres2 = (y2-mu(2)-lambda(2)*garch2)^2
series res1res2 = (y1-mu(1)-lambda(1)*garch1)*(y2-mu(2)-lambda(2)*garch2)
' LOG LIKELIHOOD
logl bvgarch
'old values
'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))
' squared errors and cross errors
bvgarch.append @logl logl
bvgarch.append sqres1 = (y1-mu(1)-lambda(1)*garch1)^2
bvgarch.append sqres2 = (y2-mu(2)-lambda(2)*garch2)^2
bvgarch.append res1res2 = (y1-mu(1)-lambda(1)*garch1)*(y2-mu(2)-lambda(2)*garch2)
' calculate the variance and covariance series
bvgarch.append var_y1 = omega(1)^2 + beta(1)^2*var_y1(-1) + alpha(1)^2*sqres1(-1)
bvgarch.append var_y2 = omega(3)^2 + omega(2)^2 + beta(2)^2*var_y2(-1) + alpha(2)^2*sqres2(-1)
bvgarch.append cov_y1y2 = omega(1)*omega(2) + beta(2)*beta(1)*cov_y1y2(-1) +alpha(2)*alpha(1)*res1res2(-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, m=100, c=1e-5)