'Load workfile and measure length of series
series obscount = 1
scalar obslength = @sum(obscount)

'Specify the return series
series y1 = djiaorg
series y2 = treas2

'Specify the number of iterations in the MLE (Engle & Sheppard (2001) used just one iteration)
!itermle = 1000

'Specify the initial values of the parameters
coef(1) alpha
coef(1) beta
coef(1) gamma

'Values may vary depending on which will result in the highest Likelihood value
alpha(1) = 0.04
beta(1) = 0.94
gamma(1) = 0.02

'Setting the sample
sample s0    @first+1 @last
sample s1    @first+2 @last
sample sf    @first+3 @last
sample sf_alt   @first+13 @last

'Initialization at sample s0
smpl s0

'Each return series is modeled with their respective GARCH specifications:
'Use Bollerslev-Wooldridge QML
'Standard errors

equation eq1.arch(1,1,thrsh=1,m=1000,c=1e-5,h) y1 c y1(-1) 'res_s11_1000 @ res_s11_1000
equation eq2.arch(1,1,thrsh=1,m=1000,c=1e-5,h) y2 c y2(-1) 'res_s22_1000 @ res_s22_1000

'Make residual series
eq1.makeresids e1
eq2.makeresids e2

'Make a garch series from the univariate estimates
eq1.makegarch h11
eq2.makegarch h22

'Normalizing the residuals from e to e* (named as "e1n" and "e2n")
series e1n = e1/h11^0.5
series e2n = e2/h22^0.5

'Make residual series for asymmetries in DCC model
series n1n = @recode(e1n<0,e1n*e1n,0)
series n2n = @recode(e2n<0,e2n*e2n,0)

'The Q_bar=E[en*en'] components and its sample equivalent 
series qbar11 = @mean(e1n*e1n)
series qbar12 = @mean(e1n*e2n)
series qbar21 = @mean(e2n*e1n)
series qbar22 = @mean(e2n*e2n)

'The N_bar
series nbar11 = @mean(n1n*n1n)
series nbar12 = @mean(n1n*n2n)
series nbar21 = @mean(n2n*n1n)
series nbar22 = @mean(n2n*n2n)

'The (en*en') matrix components
series e11n = e1n*e1n
series e12n = e1n*e2n
series e21n = e2n*e1n
series e22n = e2n*e2n

'The (nn*nn') matrix components
series n11n = n1n*n1n
series n12n = n1n*n2n
series n21n = n2n*n1n
series n22n = n2n*n2n

'Initialize the elements of Qt for variance targeting
series q11 = @var(e1)
series q12 = @cov(e1,e2)
series q21 = @cov(e2,e1)
series q22 = @var(e2)

'***********************************************************************************************

'Declare a Loglikelihood object
logl dcc
dcc.append @logl logl

'The elements of matrix Qt
dcc.append  q11 = (qbar11 - (alpha(1))^2*qbar11 - (beta(1))^2*qbar11 + (gamma(1))^2*nbar11) + (alpha(1))^2*e11n(-1) + (gamma(1))^2*n11n(-1) + (beta(1))^2*q11(-1)
dcc.append  q12 = (qbar12 - (alpha(1))^2*qbar12 - (beta(1))^2*qbar12 + (gamma(1))^2*nbar12) + (alpha(1))^2*e12n(-1) + (gamma(1))^2*n12n(-1) + (beta(1))^2*q12(-1)
dcc.append  q21 = (qbar21 - (alpha(1))^2*qbar21 - (beta(1))^2*qbar21 + (gamma(1))^2*nbar21) + (alpha(1))^2*e21n(-1) + (gamma(1))^2*n21n(-1) + (beta(1))^2*q21(-1)
dcc.append  q22 = (qbar22 - (alpha(1))^2*qbar22 - (beta(1))^2*qbar22 + (gamma(1))^2*nbar22) + (alpha(1))^2*e22n(-1) + (gamma(1))^2*n22n(-1) + (beta(1))^2*q22(-1)

'As input to detQQQ
dcc.append  q12n = ((qbar11 - (alpha(1))^2*qbar11 - (beta(1))^2*qbar11 + (gamma(1))^2*nbar11) + (alpha(1))^2*e11n(-1) + (gamma(1))^2*n11n(-1) + (beta(1))^2*q11(-1))/((abs(q11)^0.5)*(abs(q22)^0.5))
dcc.append  q21n = ((qbar21 - (alpha(1))^2*qbar21 - (beta(1))^2*qbar21 + (gamma(1))^2*nbar21) + (alpha(1))^2*e21n(-1) + (gamma(1))^2*n21n(-1) + (beta(1))*2*q21(-1))/((abs(q22)^0.5)*(abs(q11)^0.5))

'Setting up the Loglikelihood function, assuming that resid ~ N(0,H)

'The Loglikelihood function is L' = -0.5*{summation from t=1 to T of ([log of determinant of inverse(diag[Qt])*Qt*inverse(diag[Qt])] + [en'*inverse(diag[Qt])*Qt*inverse(diag[Qt])*en])

'Taking the adjoint{inverse(diag[Qt])*Qt*inverse(diag[Qt])}
dcc.append  detQQQ = 1 - q12n*q21n

'Tomando el adjoint{inverse(diag[Qt])*Qt*inverse(diag[Qt])}
dcc.append  cofact11 =   1*1
dcc.append  cofact12 =(-1)*q21n
dcc.append  cofact21 =(-1)*q12n
dcc.append  cofact22 =   1*1

'Taking the inverse{inverse(diag[Qt])*Qt*inverse(diag[Qt])}
dcc.append  invQQQ11 = cofact11/detQQQ
dcc.append  invQQQ12 = cofact12/detQQQ
dcc.append  invQQQ21 = cofact21/detQQQ
dcc.append  invQQQ22 = cofact22/detQQQ

'Taking the en'*inverse{inverse(diag[Qt])*Qt*inverse(diag[Qt])}*en
dcc.append  enQQQen11 = e1n*invQQQ11*e1n
dcc.append  enQQQen12 = e1n*invQQQ12*e2n
dcc.append  enQQQen21 = e2n*invQQQ21*e1n
dcc.append  enQQQen22 = e2n*invQQQ22*e2n

'Append the loglikelihood function
'Instead of log(detQQQ) use log(abs(detQQQ))
dcc.append logl = -0.5*(log(abs(detQQQ)) + (enQQQen11+enQQQen21+ enQQQen12+enQQQen22))

'Specifies the sample data where the estimation will be made
smpl sf

'Estimates the parameters now using BHHH algorithm
dcc.ml(b, showopts, m=!itermle, c=1e-5) 

'A detQQQnpd = 0 indicates detQQQ is positive definite
series count = (detQQQ<=0)
scalar detQQQnpd = @sum(count)

'Display the estimated parameters
show dcc.output

'Forecast the q's by initializing them
series  q11f = 0
series  q12f = 0
series  q21f = 0
series  q22f = 0

'Specify the sample period to be forecasted
smpl sf
series  q11f = (qbar11 - (alpha(1))^2*qbar11 - (beta(1))^2*qbar11 + (gamma(1))^2*nbar11) + (alpha(1))^2*e11n(-1) + (gamma(1))^2*n11n(-1) + (beta(1))^2*q11(-1)
series  q12f = (qbar12 - (alpha(1))^2*qbar12 - (beta(1))^2*qbar12 + (gamma(1))^2*nbar12) + (alpha(1))^2*e12n(-1) + (gamma(1))^2*n12n(-1) + (beta(1))^2*q12(-1)
series  q21f = (qbar21 - (alpha(1))^2*qbar21 - (beta(1))^2*qbar21 + (gamma(1))^2*nbar21) + (alpha(1))^2*e21n(-1) + (gamma(1))^2*n21n(-1) + (beta(1))*2*q21(-1)
series  q22f = (qbar22 - (alpha(1))^2*qbar22 - (beta(1))^2*qbar22 + (gamma(1))^2*nbar22) + (alpha(1))^2*e22n(-1) + (gamma(1))^2*n22n(-1) + (beta(1))^2*q22(-1)

'To plot the time-varying conditional correlation
smpl sf_alt
series rt = q12f/(@sqr(q11f)*@sqr(q22f))
graph rtgraph.line rt
show rtgraph

delete e1n e2n
delete q12n q21n
delete n1n n2n
delete n11n n12n n21n n22n
delete e11n e12n e21n e22n
delete qbar*
delete xbar*
delete cofact*
delete enqqqen*
delete invqqq*
delete enqqqen*
delete eigenvect*
delete invqqq*

show detqqqnpd

series eigmins
scalar eigmin

for !j = 4 to obslength
   'Input the hijsim to the elements of a 2x2 Q matrix
   sym(2) Q
   Q(1,1) = q11f(!j)
   Q(2,1) = q21f(!j)
   Q(2,2) = q22f(!j)

   'Take the eigenvalues of Q
   vector(4) Qeigenval
   Qeigenval = @eigenvalues(Q)

   'Collect the min eigenvalues in each Qt
   eigmins(!j) = @min(Qeigenval)
next

eigmin = @min(eigmins)
show eigmin