* =====================================
* Treshold-GARCH Modelling
* Ing. Tom Heryán, Ph.D.
* =====================================


**# Bookmark #1
* =====================================
// I. // P Ř E P S A T   Z D R O J 
clear
import excel "K:\INS_Tommík_2023.xlsx", sheet("DATA_all") firstrow
* =====================================


// elegant is business calendar - now that bug with tin( ) has been fixed
bcal create buscal, from(Date) replace
bcal load buscal
generate bcal_Date = bofd("buscal", Date)
assert !missing(bcal_Date) if !missing(Date)
format %tbbuscal bcal_Date
tsset bcal_Date











* =====================================
// II. // U P R A V I T   P R O M Ě N N É
* =====================================

encode dummy, generate(dummy2)		// string na kategorické proměnné
destring DAX_Open, replace force float
destring DAX_Close, replace force float

tsset Date, daily			// TOTO VELMI DŮLEŽITÉ !!!

label variable DAX_Close "DAX (prává osa)"
label variable BMW_Close "BMW"
label variable VW_Close "Volkswagen"
label variable RWE_Close "RWE"

// 2 // Denní výnos na akcii, grafická analýza, stacionarita

generate rDAX = (log(DAX_Close)-log(DAX_Open))*100
generate rBMW = (log(BMW_Close)-log(BMW_Open))*100
generate rVW = (log(VW_Close)-log(VW_Open))*100
generate rRWE = (log(RWE_Close)-log(RWE_Open))*100

label variable rDAX "DAX"
label variable rBMW "BMW"
label variable rVW "Volkswagen"
label variable rRWE "RWE"










**# Bookmark #2
* =====================================
// III. // NÍŽE  N A H R A D I T  DVĚ KLÍČOVÁ SLOVA 
  
* A) " close"  za  " XXX_Close" (!! POZOR mezera před)
* B) "returns"  za  "rXXX" (individuální názvy obou výše)

* =====================================



// plot the index - reproduce graph CLOSE 
twoway (line close Date), ///
		xlabel(, format(%tdCCYY)) xtitle("") ylabel(,angle(0)) ytitle("") ///
		name(index, replace) scheme(s2mono)
 
graph export volatility1.pdf, as(pdf) replace 


// plot the percentage returns - reproduce grapg RETURNS
twoway (line returns Date), ///
		xlabel(, format(%tdCCYY)) xtitle("") ylabel(,angle(0)) ytitle("") ///
		name(rets, replace) scheme(s2mono)

graph export "volatility2.pdf", as(pdf) replace 


// plot correlograms - reproduce Figure CORRELOGRAMS
ac returns, lags(20) title("ACF returns", size(medsmall)) ylabel(,angle(0)) ///
			xtitle("Lags") ytitle("Autocorrelations") recast(bar) ///
			fcolor(none)  note("") lcolor(black) lwidth(medthick) name(ac1,replace) scheme(s2mono)			
generate returns2 = returns^2

ac returns2, lags(20) title("ACF squared returns", size(medsmall)) ylabel(,angle(0)) ///
			xtitle("Lags") ytitle("Autocorrelations") recast(bar) ///
			fcolor(none)  note("") lcolor(black) lwidth(medthick) name(ac2,replace) scheme(s2mono)
			
graph combine ac1 ac2, cols(1) imargin(tiny) scheme(s2mono)

graph export "volatility3.pdf", as(pdf) replace 

*****************************************************

/* ----------------------- GARCH(1,1) estimation ----------------------- */

// estimate a GARCH(1,1) model - reproduce Table 2
arch returns, arch(1) garch(1) vce(robust) distribution(t) vsquish nolog cformat(%6.4f)
estimates store mod1
predict double h, var
generate ha = sqrt(252*h)
predict double resids, res
generate stres = resids/sqrt(h)
generate stres2 = stres^2

corrgram stres2, lag(20)

// SIC says which GARCH(p,q) is the best (LOWEST SIC)

quietly {
   noisily display  "p" _col(12)  "q" _col(20) "SIC"
   noisily display  _dup(30) "-"
      forvalues p = 1/5 {
         forvalues q = 1/2 {
            arch returns, arch(1/`p') garch(1/`q')
            estat ic
            mat stats = r(S)
            noisily display `p' _col(12) `q'  _col(20) stats[1,6]
         }
      }
}

// arch returns, arch(1/2) garch(1/2) nolog vsquish cformat(%6.4f) 
// !!! FOR GARCH (p,q) via SIC


/* ----------------------- Mean reversion in vol  ----------------------- */


// is the sum of α and β significantly less than one
// estimates restore mod1
estimates restore mod1
estimates replay mod1
nlcom _b[ARCH:L.arch]+_b[ARCH:L.garch]
test _b[ARCH:L.arch]+_b[ARCH:L.garch] = 1

nlcom sqrt(252*(_b[ARCH:_cons]/(1-_b[ARCH:L.arch]-_b[ARCH:L.garch])))

// display "Annualized volatility    " sqrt(252*(_b[ARCH:_cons]/(1-_b[ARCH:L.arch]-_b[ARCH:L.garch])))


  

// estimating half-life
display "Estimated half-life   " ceil(log(0.5)/log(_b[ARCH:L.arch]+_b[ARCH:L.garch]))


capt drop resids
capt drop h
capt drop ha
predict double resids, res
predict double h, var
generate ha = sqrt(252*h)


// reproduce Figure 4
twoway (line ha Date), ///
		xlabel(, format(%tdCCYY)) xtitle("") ylabel(,angle(0)) ytitle("") ///
		title("Annualized conditional variance") ///
		name(vol, replace)  scheme(s2mono)
		
graph export "volatility4.pdf", as(pdf) replace 

*****************************************************

// estimate a TARCH(1,1) model - reproduce Table 4
// Engle and Patton say they use robust standard errors -- but do they????
// note sign change because of specification
// use t distribution
arch returns, arch(1) garch(1) tarch(1) vce(robust) distribution(t) nolog vsquish cformat(%6.4f)

// estimate GARCH(1,1)-X model === ADD VARIABLE INTO VARIANCE EQ.
// note Stata specfication is different to EVIEWS het() is exponentiated
// arch djrets, arch(1) garch(1) het(EURIBOR) vce(robust) vsquish nolog cformat(%6.4f)

*****************************************************

estimates clear
// estimate models
arch returns, arch(1) garch(1) vce(robust) distribution(t) nolog
estimates store mod1
arch returns, arch(1) garch(1) tarch(1) vce(robust) distribution(t) nolog
estimates store mod2
// arch returns, arch(1) garch(1) het(EURIBOR) vce(robust) distribution(t) nolog
// estimates store mod3

estimates table mod*, b(%6.4f) se(%6.4f) stats(ll) stfmt(%6.2f)

END