I am trying to program the estimation of a GARCH-in-mean model with DCC. The code I have is for a bivariate (asymmetric) GARCH with DCC and I do not know how to program the feedback into the mean equation. Any help would be greatly appreciated.
Code: Select all
' ORIGINAL
'H = (Dt)*(Rt)*(Dt)
'Rt = inverse(diag[Qt])*Qt*inverse(diag[Qt])
'Qt = (Qbar - alpha*Qbar - beta*Qbar) + alpha*(en*en') + beta*Qt
'This is a DCC with 2 series, a DCC(1,1) Model
' [0] Load workfile and measure length of series
load g:\siri\dcc.wf1
series obscount = 1
scalar obslength = @sum(obscount)
' [1] SPECIFY THE ASSET RETURN SERIES
'For two return series
series y1 = ret_totmt
series y2 = ret_5x
' [2] SPECIFY THE NUMBER OF ITERATIONS IN THE MLE, Engle and Sheppard (2001) used just one iteration
!itermle = 100
' [3] SPECIFY THE INITIAL VALUES OF THE PARAMETERS
'Coefficient vectors
coef(1) alpha
coef(1) beta
'VARY these depending on which will result in the highest likelihood value
alpha(1) = 0.011390
beta(1) = 0.982627
'Setting the sample
sample s0 @first+1 @last
sample s1 @first+2 @last
sample sf @first+3 @last
sample sf_alt @first+13 @last
'sample s0 @first+561 @last-590
'sample s1 @first+562 @last-590
'sample sf @first+563 @last-590
'sample sf_alt @first+573 @last-590
'Initialization at sample s0
smpl s0
'Each return series is modeled with their respective mean equations and GARCH specifications:
'Use Bollerslev-Wooldridge QML
'standard errors
equation eq1.arch(egarch,m=100,c=1e-5,h) y1 c res_b11_100 @ res_b11_100
equation eq2.arch(thrsh=1,m=100,c=1e-5,h) y2 c y2(-1) res_s_1000 @ res_s_1000
'Make residual series
eq1.makeresids e1
eq2.makeresids e2
'Make a garch series from the univariate estimates
'p.202, p.336 in Command Reference
eq1.makegarch h11
eq2.makegarch h22
'Normalized the resid e, e* -> en, wrt to h^0.5
series e1n = e1/h11^0.5
series e2n = e2/h22^0.5
'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 (en*en') matrix components
series e11n = e1n*e1n
series e12n = e1n*e2n
series e21n = e2n*e1n
series e22n = e2n*e2n
'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)*qbar11 - beta(1)*qbar11 + alpha(1)*e11n(-1) + beta(1)*q11(-1)
dcc.append q12 = qbar12 - alpha(1)*qbar12 - beta(1)*qbar12 + alpha(1)*e12n(-1) + beta(1)*q12(-1)
dcc.append q21 = qbar21 - alpha(1)*qbar21 - beta(1)*qbar21 + alpha(1)*e21n(-1) + beta(1)*q21(-1)
dcc.append q22 = qbar22 - alpha(1)*qbar22 - beta(1)*qbar22 + alpha(1)*e22n(-1) + beta(1)*q22(-1)
'As input to detQQQ
dcc.append q12n = (qbar12 - alpha(1)*qbar12 - beta(1)*qbar12 + alpha(1)*e12n(-1) + beta(1)*q12(-1))/((abs(q11)^0.5)*(abs(q22)^0.5))
dcc.append q21n = (qbar21 - alpha(1)*qbar21 - beta(1)*qbar21 + alpha(1)*e21n(-1) + beta(1)*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])
'So the determinant of inverse(diag[Qt])*Qt*inverse(diag[Qt])
dcc.append detQQQ = 1 - q12n*q21n
'Taking the 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 s1
'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)*qbar11 - beta(1)*qbar11 + alpha(1)*e11n(-1) + beta(1)*q11(-1)
series q12f = qbar12 - alpha(1)*qbar12 - beta(1)*qbar12 + alpha(1)*e12n(-1) + beta(1)*q12(-1)
series q21f = qbar21 - alpha(1)*qbar21 - beta(1)*qbar21 + alpha(1)*e21n(-1) + beta(1)*q21(-1)
series q22f = qbar22 - alpha(1)*qbar22 - beta(1)*qbar22 + alpha(1)*e22n(-1) + beta(1)*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 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