*
* ARTIGO RBFim
*
* Replication file for Chan and Maheu(2002), "Conditional Jump Dynamics
* in Stock Market Returns", JBES, vol 20, no. 3, 377-389.
*
* ARJI model (table 3)
*
source C:\Users\Max\Desktop\Artigos\SBFin_2019\PETROLEO\jumpgarch.src
*
open data C:\Users\Max\Desktop\Artigos\SBFin_2019\PETROLEO\spot_2.txt
data(format=free,org=columns) 1 2008 year month day close logret
cal(julian=%julianfromymd(year,month,day))
*
set r = 100.0*logret
stats r
*
sstats(max) / %if(year(t)<=2017,t,0)>>x2017
sstats(min) / %if(year(t)>=2010,t,1000000)>>x2010
compute end2017=fix(x2017),start2010=fix(x2010)
*
set crisis = t>=2010:4:10.and.t<=2010:5:10.or.$
             t>=2011:1:10.and.t<=2011:8:1.or.$
             t>=2012:5:1.and.t<=2012:6:30.or.$
             t>=2014:10:10.and.t<=2015:6:1
*
*
graph(shading=crisis, STYLE=LINE) 1
# CLOSE
graph(shading=crisis, STYLE=LINE) 1
# r
*
*
spgraph(hfields=2,footer="Figura 0 ")
graph(shading=crisis,header="(a) Cotações US$/barril",style=lines)
# close
graph(shading=crisis,header="(b) Retornos (%)",style=lines)
# r
spgraph(done)

* ARCH Test
linreg r / resids
# constant
@archtest(lags=6,form=lm,span=1) resids
*
* Test autocorrelation Ljung-Box Q
correlate(number=12,qstats) r
*
*
set u   = 0.0
set h   = 0.0
*
* We need these series to form the recursive estimates of the jump intensity
*
set lambda_t  = 0.0
set xi_t      = 0.0
*
* These save the (possibly) time-varying mean and variance of the jump process.
*
set theta_t   = 0.0
set deltasq_t = 0.0
*
* Define formula for the mean model and set up its parameter set
*
nonlin(parmset=meanparms) mu a1
frml uf = r-mu-a1*r{1}
*
* Define formula for the base GARCH variance and set up its parameter set
*
nonlin(parmset=garchparms) omega alpha beta
frml hf = omega+alpha*u{1}^2+beta*h{1}
*
* Define the ARJI formula for the intensity.
*
declare real lambda0 rho gamma
frml lambdaf = lambda0+rho*lambda_t{1}+gamma*xi_t{1}
****************************************************************************
*
* The parmsets and logl definitions need to change from model to model
*
****************************************************************************
*
* Constant intensity jump model through 2017 (column 1 in table 3)
*
garch(p=1,q=1,reg,resids=u,hseries=h) 4 end2017 r
# constant r{1}
compute ssvar=%beta(%nregmean+1)/(1-%beta(%nregmean+2)-%beta(%nregmean+3))
set u * %regstart()-1 = 0.0
set h * %regstart()-1 = ssvar
compute mu=%beta(1),a1=%beta(2),a2=%beta(3)
compute omega=%beta(%nregmean+1),alpha=%beta(%nregmean+2),beta=%beta(%nregmean+3)
nonlin(parmset=poissonparms) lambda0
nonlin(parmset=jumpparms) zeta0 eta0
compute eta0=0.0,zeta0=sqrt(ssvar),lambda0=.1
frml logl = u=uf,h=hf,lambda_t=lambda0,deltasq_t=zeta0^2,theta_t=eta0,$
  arjigarch(u,h,lambda_t,deltasq_t,theta_t,xi_t)
maximize(parmset=meanparms+garchparms+jumpparms+poissonparms,$
 pmethod=simplex,piters=2,method=bfgs,$
 title="Constant Intensity Jump GARCH, Full Sample") logl 4 end2017
*
* TESTE DE DIAGNOSTICO PARA OQUADRADO DOS RESIDUOS
*
* Diagnostics on squared standardized residuals (taking into account the
* mean and variance of the Poisson process) and on the jump intensity
* residuals.
*
*
set ustdsq = (u-lambda_t*theta_t)^2/(h+lambda_t*(theta_t^2+deltasq_t))
@westchotest(title="Q Test on Standardized Squared Residuals",number=15) ustdsq %regstart() %regend()
@westchotest(title="Q Test on Jump Intensity Residuals",number=15) xi_t %regstart() %regend()
*
* Save this for later use
*
compute lambdafixed=lambda0
****************************************************************************
*
* ARJI - GARCH
*
* Simple ARJI model (column 2 in table 3). We'll just take the converged
* estimates from the constant model as our initial guesses.
*
nonlin(parmset=poissonparms) lambda0 rho gamma rho=gamma
compute lambda0=.021,rho=.25,gamma=0.1
frml logl = u=uf,h=hf,lambda_t=lambdaf,deltasq_t=zeta0^2,theta_t=eta0,$
   ARJIgarch(u,h,lambda_t,deltasq_t,theta_t,xi_t)
maximize(parmset=meanparms+garchparms+jumpparms+poissonparms,$
 pmethod=simplex,piters=2,method=bfgs,iterations=50,title="ARJI-GARCH") logl 4 end2017
*
*
* TESTE DE DIAGNOSTICO PARA OQUADRADO DOS RESIDUOS
*
set ustdsq = (u-lambda_t*theta_t)^2/(h+lambda_t*(theta_t^2+deltasq_t))
@westchotest(title="Q Test on Standardized Squared Residuals",number=15) ustdsq %regstart() %regend()
@westchotest(title="Q Test on Jump Intensity Residuals",number=15) xi_t %regstart() %regend()
*
*  FIGURA 1
*
set lambda_simple = lambda_t
*
set crisis = t>=2010:4:10.and.t<=2010:5:10.or.$
             t>=2011:1:10.and.t<=2011:8:1.or.$
             t>=2012:5:1.and.t<=2012:6:30.or.$
             t>=2014:10:10.and.t<=2015:6:1
*
spgraph(vfields=2,footer="Figure 1")
graph(shading=crisis,header="(a) Retornos do petróleo",style=lines)
# r
graph(shading=crisis,header="(b) Intensidade do salto",style=lines)
# lambda_simple
spgraph(done)
*
* FIGURE 2
* Figure out the entries for the particular y:m:d combinations
*
sstats(max) / $
  %if(year(t)==2012.and.month(t)==6.and.day(t)==18,t,0)>>x20120618 $
  %if(year(t)==2010.and.month(t)==4.and.day(t)==20,t,0)>>x20100420
*
* Pull out the lambda's for those entries, and compute the Poisson
* probabilities for the fixed model and the ARJI model entries.
*
compute lambda2012=lambda_simple(fix(x20120618))
compute lambda2010=lambda_simple(fix(x20100420))
set probfixed 1 11 = %poissonk(lambdafixed,t-1)
set prob2012  1 11 = %poissonk(lambda2012,t-1)
set prob2010  1 11 = %poissonk(lambda2010,t-1)
graph(footer="Figure 2. Probabilities of Number of Jumps",number=0,style=bar,key=below,$
  klabels=||"Constante","Rebaixamento de títulos","Acidente Plataforma BP"||) 3
# probfixed
# prob2012
# prob2010
*
* FIGURE 3
*
sstats(max) / $
  %if(year(t)==2012.and.month(t)==6.and.day(t)==1,t,0)>>x20120601 $
  %if(year(t)==2012.and.month(t)==6.and.day(t)==29,t,0)>>x20120629
spgraph(vfields=2,footer="Figure 3 ARJI, In-Sample R and Lambda")
scatter(style=line,header="(a) Retornos do petróleo")
# day r fix(x20120601) fix(x20120629)
scatter(style=line,header="(b) Intensidade do salto")
# day lambda_simple fix(x20120601) fix(x20120629)
spgraph(done)
****************************************************************************
*
* ARJI-R{1}^2 (column 3 in table 3)
*
nonlin(parmset=jumpparms) zeta0 zeta1 eta0 eta1 eta2
compute eta1=eta2=0.0
sstats(mean) * end2017 r^2>>meanrsq
compute zeta0=.2*zeta0,zeta1=zeta0^2/meanrsq
frml logl = u=uf,h=hf,lambda_t=lambdaf,deltasq_t=zeta0^2+zeta1*r{1}^2,$
   theta_t=eta0+eta1*%max(r{1},0.0)+eta2*%min(r{1},0.0),$
   ARJIgarch(u,h,lambda_t,deltasq_t,theta_t,xi_t)
maximize(parmset=meanparms+garchparms+jumpparms+poissonparms,$
 pmethod=simplex,piters=2,method=bfgs,title="ARJI-R{1}^2 GARCH") logl 4 end2017
*
*
* Teste de diagnóstico
*
set ustdsq = (u-lambda_t*theta_t)^2/(h+lambda_t*(theta_t^2+deltasq_t))
@westchotest(title="Q Test on Standardized Squared Residuals",number=15) ustdsq %regstart() %regend()
@westchotest(title="Q Test on Jump Intensity Residuals",number=15) xi_t %regstart() %regend()
*
*
* FIGURE 4
*
set lambda_r2 = lambda_t
*
spgraph(vfields=2,footer="Figure 4: ARJI-R{1}^2 ")
graph(shading=crisis,header="(a) Retornos dopetróleo",style=lines)
# r
graph(shading=crisis,header="(b) Intensidade do salto",style=lines)
# lambda_r2
spgraph(done)
*
* FIGURE 5
*
* Figure out the entries for the particular y:m:d
*
* sstats(max) / $
*  %if(year(t)==2012.and.month(t)==6.and.day(t)==18,t,0)>>x20120618 $
*  %if(year(t)==2010.and.month(t)==4.and.day(t)==20,t,0)>>x20100420
*
* Pull out the lambda's for those entries, and compute the Poisson
* probabilities for the fixed model and the ARJI-R{1}^2  model entries.
*
compute lambda2012=lambda_r2(fix(x20120618))
compute lambda2010=lambda_r2(fix(x20100420))
set probfixed 1 11 = %poissonk(lambdafixed,t-1)
set prob2012r2  1 11 = %poissonk(lambda2012,t-1)
set prob2010r2  1 11 = %poissonk(lambda2010,t-1)
graph(footer="Figure 5. Probabilities of Number of Jumps",number=0,style=bar,key=below,$
 klabels=||"Constante","Rebaixamento de títulos","Acidente Plataforma BP"||) 3
# probfixed
# prob2012r2
# prob2010r2
*
* FIGURE 6
*
sstats(max) / $
  %if(year(t)==2012.and.month(t)==6.and.day(t)==1,t,0)>>x20120601 $
  %if(year(t)==2012.and.month(t)==6.and.day(t)==29,t,0)>>x20120629
spgraph(vfields=2,footer="Figure 6: ARJI-R^2 and Lambda")
scatter(style=line,header="(a) Retornos do petróleo")
# day r fix(x20120601) fix(x20120629)
scatter(style=line,header="(b) Intensidade do salto")
# day lambda_r2 fix(x20120601) fix(x20120629)
spgraph(done)
*
*
*****************************************************************************
*
*
* ARJI-h(t) (column 4 in table 3)
*
nonlin(parmset=jumpparms) zeta0 zeta1 eta0 eta1 eta2
frml logl = u=uf,h=hf,lambda_t=lambdaf,deltasq_t=zeta0^2+zeta1*h,$
  theta_t=eta0+eta1*%max(r{1},0.0)+eta2*%min(r{1},0.0),$
  ARJIgarch(u,h,lambda_t,deltasq_t,theta_t,xi_t)
maximize(parmset=meanparms+garchparms+jumpparms+poissonparms,$
 pmethod=simplex,piters=2,method=bfgs,title="ARJI-h GARCH") logl 4 end2017
*
*
*
set ustdsq = (u-lambda_t*theta_t)^2/(h+lambda_t*(theta_t^2+deltasq_t))
@westchotest(title="Q Test on Standardized Squared Residuals",number=15) ustdsq %regstart() %regend()
@westchotest(title="Q Test on Jump Intensity Residuals",number=15) xi_t %regstart() %regend()
*
* FIGURE 7
*
set lambda_h = lambda_t
*
spgraph(vfields=2,footer="Figure 7: ARJI-h(t) ")
graph(shading=crisis,header="(a) Retornos do petróleo",style=lines)
# r
graph(shading=crisis,header="(b) Intensidade do Salto",style=lines)
# lambda_h
spgraph(done)
*
* FIGURE 8
* Figure out the entries for the particular y:m:d
*
sstats(max) / $
  %if(year(t)==2012.and.month(t)==6.and.day(t)==18,t,0)>>x20120618 $
  %if(year(t)==2010.and.month(t)==4.and.day(t)==20,t,0)>>x20100420
*
* Pull out the lambda's for those entries, and compute the Poisson
* probabilities for the fixed model and the ARJI model entries.
*
compute lambda2012=lambda_h(fix(x20120618))
compute lambda2010=lambda_h(fix(x20100420))
set probfixed 1 11 = %poissonk(lambdafixed,t-1)
set prob2012h  1 11 = %poissonk(lambda2012,t-1)
set prob2010h  1 11 = %poissonk(lambda2010,t-1)
graph(footer="Gráfico 1. Probabilities of Number of Jumps",number=0,style=bar,key=below,$
  klabels=||"Constante","Rebaixamento de títulos","Acidente Plataforma BP"||) 3
# probfixed
# prob2012h
# prob2010h
*
*
* FIGURE 9
*
sstats(max) / $
  %if(year(t)==2012.and.month(t)==6.and.day(t)==1,t,0)>>x20120601 $
  %if(year(t)==2012.and.month(t)==6.and.day(t)==29,t,0)>>x20120629
spgraph(vfields=2,footer="Figura 9: ARJI-h(t) and Lambda")
scatter(style=line,header="(a) Retornos do petróleo")
# day r fix(x20120601) fix(x20120629)
scatter(style=line,header="(b) Intensidade do salto")
# day lambda_h fix(x20120601) fix(x20120629)
spgraph(done)
*
*