# Risk Management under Volatility Clustering

## Introduction

Here are a few things covered during this week.

• What is the objective of risk management?
• To create a statistical model for future returns , in particular the left tail.
• Model if often based on historical returns ( the only thing that we can observe)
• Future Distribution of Returns(TR) = Function (Historical Returns , Some assumptions)
• Modelling assumptions so far
• Future Distribution ( WILSHIRE 5000/ index fund) == Historical Distribution ( True for some time but not always)
• Currently our distribution of choice is re-scaled T-distribution as it models the data better than normal distribution .
• We have also simulated observed distribution from the data
• Parameters of historical distribution ( mean, sigma and/ or df) are calculated without paying attention to the ordering of the data
• Implies, even if we randomly permeates logret series distribution parameter will still be same.
• What if we have some information in the past which can inform us about the future distribution for logret? In other words, what is the impact of serial correlation if any.
• Serial correlation
• Positive serial correlation - an above average return ( > mean) tend to be followed by another above - average return
• Negative serial correlation - above average return is followed by below average return
• No correlation / random wal k - an above average return does not increase the likelihood of being followed by an above avg or below avg return.
• There was a long history of research in financial markets that look at serial correlation. Now to understand what serial correlation means,
• Correlation - let’s start with the correlation between two variables, let’s say, X and Y. The correlation coefficient between X and Y measures how they co-move together. We say that X and Y are positively correlated when they tend to move up and down together. To be more precise, if X is above its mean, then Y is likely to be also above its mean.
• Serial correlation - How a variable X and its own past are correlated?
• Significance of serial correlation in Finance is associated with the concept of Market Efficiency.
• The theory of market efficiency says that the observed price of an asset such as a stock, fully reflects all information about that asset

• When new information comes in, it will change the price of the asset. Good news presumably will lead to an increase in its price. Conversely, bad news will lead to a decrease in its price.

• Now, because all existing information has been incorporated in the price of an asset, any new information must be a surprise.

• If the surprise is good news, the price of the asset will go up, and if the surprise is bad news, the price of the asset will go down. But there is no way to tell ahead of time whether the news is good or bad. So there is no way to tell if the price will go up or down.

• This is so called random walk model of prices.

• This random walk model just says there is no way to predict if the future price of an asset will go up or go down relative to its current price.

• An implication of the random walk model is that returns have no serial correlation. So testing for serial correlation in stock return has been viewed as testing for market efficiency

• Motivation & Definition for Volatility Clustering - In financial time series data, it is often observed
• Large returns +ve/ -ve negative tend to be followed by more large returns

• High Volatility days are followed by more High Volatility days

• Low Volatility days are followed by more Low Volatility days

• This phenomenon is called Volatility clustering

• It is calculated mathematically by measuring autocorrelation of absolute values of log returns

• A time series exhibiting Volatility clustering implies ordering of the data is important

• Strong volatility clustering implies we can predict volatility of stock prices for future values based on the past values.

• This is exploited in GARCH model ( Generalized Auto-regressive Conditional Heteroskedasticity)

• We can test volatility clustering/ effect of ordering by

• Plotting acf on absolute values of logreturns - ( Should show strong evidence of serial correlation )

• Shuffling/ Permeating returns and apply same logic again ( Should show no evidence of serial correlation on shuffled series)

• GARCH model
• WILSHIRE 5000 ACF plot indicates no/ weak evidence of serial correlation => Difficult to predict mean of future values based on historical values
• WILSHIRE 5000 ABS(ACF) plot indicates strong volatility clustering =>Volatility / Standard Deviation is changing over time in a predictable manner.
• What we want is a statistical model which allows us to predict future volatility- GARCH
• GARCH model has many extensions. We are going to focus on GARCH(1,1). It is one of the simplest GARCH model .But it does a very good job on WILSHIRE 5000 data in modelling volatility clustering.
• GARCH(1,1)-normal model is equivalent to normal model when we assume constant variance.
• Equations

\begin{align} & \text{Mean Equation: }r_t = a_0 + \sqrt{h_t}\epsilon_t \\ & \text{Variance Equation: } h_t = \alpha_0 + \beta_1 h_{t-1} + \alpha_1\epsilon_{t-1} \\ &\text{Distribution Equation: } \epsilon_t \sim N(0,1) \end{align}

• Notations:

• $a_0$ : Expected Returns
• $\sqrt{h_t}\epsilon_t$ : Unexpected Returns
• $r_t$ : Total returns of series with time varying volatility
• $h_t$ : predictable variance of $r_t$ changing over time
• $\epsilon_t$ : Normal Distribution
• If $h_t$ is constant above equation reduces to normal model for logret.

• For Volatility clustering to happen: $\alpha_1 > 0 \text{ and } \beta_1 >0$

• GARCH(1,1)-t model is equivalent to normal model when we assume constant variance.
• Equations

\begin{align} & \text{Mean Equation: }r_t = a_0 + \sqrt{h_t}\epsilon_t \\ & \text{Variance Equation: } h_t = \alpha_0 + \beta_1 h_{t-1} + \alpha_1\epsilon_{t-1} \\ &\text{Distribution Equation: } \epsilon_t \sim \frac{t(v)}{\sqrt{\frac{v}{v-2}}} \end{align}

• Here we replace chosen distribution (norm) with re-scaled -t model

• 1-day ahead garch model can be simulated using ugarchboot function.
• VaR and ES from garch take time ordering into account
• When Fitted values are low => VaR and ES are also low.( Only apply for 1-day)
• Why is GARCH model better?
• When volatility is high => Risk is also high => Portfolio manager can take some action to mitigate risk
• When volatility is low => Risk is low => No action as risk reduction may be costlier to portfolio.
• VaR graph can be produced by using ugarchroll function.
• Summary
• Basics of Financial risk management
• Analyze logret of portfolio
• Distribution: typically non-normal
• We also have volatility clustering : a behavior typical of asset returns
• GARCH(1,1)-t can explain both
• VaR and ES change over time
• Ideas can be applied to different types of portfolios
• Margin borrowing: $10 K stocks with •$5K your money
• $5k from a brokerage firm • This is a margin position • You are hoping to make money from stock price going up • Risk -> Price can go down. How much money can you afford to loose -> ES over a time period( like 10 days) with a certain confidence level( like 95%) • Stock/Bond Portfolio: • Retirement Account • 50% US Equities • 50% long term corporate bond • Plan to retire in 1 year • No plan to add more money or take money out next year • How much can you loose in 1 year? • Can you afford to loose this amount of money for retirement in 1 year? ## Imports library(quantmod) Loading required package: xts Loading required package: zoo  Attaching package: 'zoo' The following objects are masked from 'package:base': as.Date, as.Date.numeric Loading required package: TTR Registered S3 method overwritten by 'quantmod': method from as.zoo.data.frame zoo  library(ggplot2) library(purrr) # partial function library(glue) # string interpolation library(moments) library(MASS) library(metRology)  Attaching package: 'metRology' The following objects are masked from 'package:base': cbind, rbind library(rugarch) Loading required package: parallel  Attaching package: 'rugarch' The following object is masked from 'package:purrr': reduce The following object is masked from 'package:stats': sigma Set Seed. RNGkind(sample.kind="Rounding") set.seed(123789) ## Functions Show the code filterSeries <- function(s, date_filter="1979-12-31/2017-12-31", new_name='TR'){ s <- na.omit(s); s s <- s[date_filter]; s names(s) <- new_name; s return(s) } getSeries <- function(name, date_filter="1979-12-31/2017-12-31", src='FRED', new_name='TR', auto.assign=FALSE){ s <- getSymbols(name, src=src, auto.assign = auto.assign); s s <- filterSeries(s) return(s) } get_df <- function(tbl){ df <- data.frame(date=index(tbl), TR=as.numeric(tbl$TR))
return(df)
}

logret_1 <- function(tbl){
return(diff(log(tbl))[-1])
}

aggr_fn <-  function(tbl, fn){
result <- fn(tbl, sum)
return(result)
}

logret_w <- partial(aggr_fn, fn = apply.weekly)
logret_m <- partial(aggr_fn, fn = apply.monthly)
logret_q <- partial(aggr_fn, fn = apply.quarterly)
logret_a <- partial(aggr_fn, fn = apply.yearly)

logret_ndays <- function(tbl, days=10){
len_tbl <-  (nrow(tbl)-days+1)
new_TR <- sapply(seq(1:len_tbl), function(i){
sum(tbl[i:(i+days-1)])
})
tbl[1:len_tbl]$TR <- new_TR return(tbl[1:len_tbl]) } ret <- function(logret){ return(exp(logret)-1) } get_norm_series <- function(mu, sig, n=100000, set_seed=TRUE){ if (set_seed){ RNGkind(sample.kind="Rounding") set.seed(123789) } return(rnorm(n,mu, sig)) } get_sim_series <- function(s, n=100000, set_seed=TRUE){ if (set_seed){ RNGkind(sample.kind="Rounding") set.seed(123789) } return(sample(as.vector(s), n, replace=TRUE)) } get_t_series <- function(mu, sig, df, n=100000, set_seed=TRUE){ if (set_seed){ RNGkind(sample.kind="Rounding") set.seed(123789) } return(rt.scaled(n,mean=mu,sd=sig,df=df)) } get_norm_series_from_s <- function(s, n=100000, set_seed=TRUE){ mu = mean(s) sig = sd(s) return(get_norm_series(mu, sig, n=n, set_seed = set_seed)) } get_fit_params <- function(s){ fitpars <- s |> as.vector() |> fitdistr(densfun='t') return(fitpars) } get_t_series_from_s <- function(s, n=100000, fitpars=NULL, set_seed=TRUE){ if(is.null(fitpars)){ fitpars <- s |> get_fit_params() } t.fit <- fitpars return(get_t_series(mu=t.fit$estimate,
sig=t.fit$estimate, df=t.fit$estimate,
n=n,
set_seed=set_seed))
}

get_sim_VaR <- function(s, alpha=0.05){
return(quantile(s, alpha) |> unname())
}

get_sim_ES <- function(s, VaR){
# VaR  <- get_sim_VaR(s, alpha=alpha)
return(mean(s[s<VaR]))
}

get_VaR <- function(tbl,alpha=0.05){
mu = mean(tbl)
sig = sd(tbl)
VaR = qnorm(alpha, mean=mu, sd=sig) #quantile of normal distribution(only not for any other distribution)
return(VaR)
}

get_ES <- function(tbl, alpha=0.05){
mu = mean(tbl)
sig = sd(tbl)
ES = mu - sig*dnorm(qnorm(alpha, mean=0, sd=1), 0, 1)/alpha #quantile of normal distribution(only not for any other distribution)
return(ES)
}

expected_loss <- function(portfolio, risk_num){ #95% confidence that we won't loose more than this amount
return(portfolio*(exp(risk_num)-1))
}

err <- function(a, b){
return((a-b)/a*100)
}

get_dist_plot <- function(tbl){
mu = mean(tbl)
sig = sd(tbl)
p <- get_df(tbl) |>
ggplot(mapping = aes(x=TR)) +
geom_histogram(aes(y=after_stat(density)), colour="black", fill="white",
bins=100)+
geom_density(alpha=.2, fill="violet") +
geom_vline(aes(xintercept=mu),
color="brown", linetype="dashed", linewidth=1) +
geom_vline(aes(xintercept=mu+2*sig),
color="brown", linetype="dashed", linewidth=0.1) +
geom_vline(aes(xintercept=mu-2*sig),
color="brown", linetype="dashed", linewidth=0.1) +
labs(title=glue("Distribution of TR: mu={round(mu, 8)}, sigma={round(sig, 8)}"), y='count')
return(p)
}

get_VaR_plot <- function(tbl, alpha=0.05){
mu = mean(tbl)
sig = sd(tbl)
VaR = qnorm(alpha, mean=mu, sd=sig)
p <- get_df(tbl) |>
ggplot(mapping = aes(x=TR)) +
geom_histogram(aes(y=after_stat(density)), colour="black", fill="white",
bins=100)+
geom_density(alpha=.2, fill="violet") +
geom_vline(aes(xintercept=mu),
color="brown", linetype="dashed", linewidth=1) +
geom_vline(aes(xintercept=VaR),
color="brown", linewidth=1) +
labs(title=glue("Distribution of TR: mu: {round(mu, 8)}, sigma: {round(sig, 6)}, VaR({alpha}): {round(VaR, 6)}"), y='count')
return(p)
}

get_comp_plot <- function(s, alpha=0.05){
mu <- mean(s)
sig <- sd(s)
norm_s <- get_norm_series(mu, sig)
VaR <- norm_s |> get_sim_VaR(); VaR
ES <- norm_s |> get_sim_ES(VaR=VaR); ES
s_sim <- s |> get_sim_series()
VaR_sim <-  s_sim|>  get_sim_VaR(); VaR_sim
ES_sim <- s_sim |> get_sim_ES(VaR=VaR_sim); ES_sim
err_VaR <- err(VaR, VaR_sim) |> round(digits = 2)
err_ES <- err(ES, ES_sim) |> round(digits = 2)

p <- data.frame(TR=s_sim, TR2=norm_s) |>
ggplot() +
geom_density(mapping = aes(x=TR), alpha=.2, fill="violet") +
geom_density(mapping = aes(x=TR2), alpha=.2, fill="white") +
geom_vline(aes(xintercept=ES_sim),
color="brown", linetype="dashed", linewidth=1) +
geom_vline(aes(xintercept=ES),
color="blue", linetype="dashed", linewidth=1) +
geom_vline(aes(xintercept=VaR_sim),
color="brown", linewidth=1) +
geom_vline(aes(xintercept=VaR),
color="blue", linewidth=1) +
labs(title=glue("Distribution of TR alpha({alpha}): \n
VaR_norm(blue): {round(VaR, 6)},VaR_sim(brown): {round(VaR_sim, 6)}, Error: {round(err(VaR, VaR_sim), 2)}% \n
ES_norm: {round(ES, 6)},ES_sim: {round(ES_sim, 6)}, Error: {round(err(ES, ES_sim), 2)}%"), y='count')
return(p)
}

do_analysis <-  function(s, alpha=0.05){

modelname1 <-  "Norm"
VaR1 <- s |> get_VaR() |> as.double()
ES1 <- s |> get_ES() |> as.double()
# df[nrow(df) + 1,] <- c(modelname, VaR, ES)

modelname2 <-  "NormSim"
norm_series <- s |> get_norm_series_from_s()
VaR2 <- norm_series|> get_sim_VaR(alpha=alpha) |> as.double()
ES2 <- norm_series |> get_sim_ES(VaR=VaR2) |> as.double()
# df[nrow(df) + 1,] <- c(modelname, VaR, ES)

modelname3 <-  "T-Sim"
t_series <- s |> get_t_series_from_s()
VaR3 <- t_series|> get_sim_VaR(alpha=alpha) |> as.double()
ES3 <- t_series |> get_sim_ES(VaR=VaR3) |> as.double()
# df[nrow(df) + 1,] <- c(modelname, VaR, ES)

modelname4 <-  "Sim"
sim_series <- s |> get_sim_series()
VaR4 <- sim_series|> get_sim_VaR(alpha=alpha) |> as.double()
ES4 <- sim_series |> get_sim_ES(VaR=VaR4) |> as.double()

model_names <- c(modelname1, modelname2, modelname3, modelname4)
VaRs <- c(VaR1, VaR2, VaR3, VaR4)
ESs <- c(ES1, ES2, ES3, ES4)

return(data.frame(ModelName=model_names, VaR=VaRs, ES=ESs))

}

get_multi_day_series <- function(s,
days=10,
sampling_func=get_sim_series,
n=100000,
set_seed=TRUE){

if (set_seed){
RNGkind(sample.kind="Rounding")
set.seed(123789)
}

rvec <- rep(0, n)
for( i in 1:days){
rvec <- rvec + sampling_func(s=s, n=n, set_seed=FALSE)
}
return(rvec)
}

get_multi_day_series_block <- function(s,
days=10,
n=100000,
set_seed=TRUE){

if (set_seed){
RNGkind(sample.kind="Rounding")
set.seed(123789)
}
rdat <- s |> as.vector()
rvec <- rep(0, n)
posn <- seq(from=1, to=(length(rdat)-days+1), by=1)
rpos <- sample(posn, size=n, replace=TRUE)
for( i in 1:days){
rvec <- rvec + rdat[rpos]
rpos <- rpos+1
}
return(rvec)
}

do_multi_analysis <-  function(rvec, days=10,  alpha=0.05){

modelname1 <-  glue("student-t({days})")
fitspars <- rvec |> get_fit_params() #; fitspars
sampling_func_tdist <-  partial(get_t_series_from_s,
fitpars = fitspars)
s1 <-  rvec |> get_multi_day_series(
sampling_func=sampling_func_tdist)
VaR1 <- s1 |> get_sim_VaR(alpha=alpha) |> as.double()
ES1 <- s1 |> get_sim_ES(VaR=VaR1) |> as.double()

modelname2 <-  glue("IID({days})")
s2 <-  rvec |> get_multi_day_series(days=days)
VaR2 <- s2 |> get_sim_VaR(alpha=alpha) |> as.double()
ES2 <- s2 |> get_sim_ES(VaR=VaR2) |> as.double()

modelname3 <-  glue("consecutive({days})")
s3 <-  rvec |> logret_ndays(days=days) |> get_sim_series(set_seed = TRUE)
VaR3 <- s3 |> get_sim_VaR(alpha=alpha) |> as.double()
ES3 <- s3 |> get_sim_ES(VaR=VaR3) |> as.double()

model_names <- c(modelname1, modelname2, modelname3)
VaRs <- c(VaR1, VaR2, VaR3)
ESs <- c(ES1, ES2, ES3)

return(data.frame(ModelName=model_names, VaR=VaRs, ES=ESs))

}

expected_loss_df <- function(df, portfolio=1000){
expected_loss_port = partial(expected_loss, portfolio=portfolio)
df$VaR <- sapply(df$VaR, expected_loss_port)
df$ES <- sapply(df$ES, expected_loss_port)
return(df)
}

## Future Vs Historical Distribution

wilsh <- getSeries("WILL5000IND")
wilsh |> head(2)
             TR
1979-12-31 1.90
1980-01-02 1.86
wilsh |> tail(2)
               TR
2017-12-28 124.33
2017-12-29 123.67
wilsh |> logret_1() |> acf() ACF Intrepretation
• Dashed lines represent 95% confidence interval for ACF

• If many vertical lines are outside dashed band => We have evidence of serial correlation.

• If all or most lines are inside dashed bands => no strong evidence of serial correlation.

• Above plot shows , We do not find much evidence of serial correlation in the daily log returns of the Wilshire 5000 index. That means the direction of stock prices are not predictable.

## Volatility Clustering

wilsh |> logret_1() |> abs() |> acf() All lines for autocorrelation on absolute values are above dashed band showing strong possibility of serial correlation and thus volatility clustering

wilsh |> logret_1() |> abs() |> as.vector() |> sample(size=nrow(wilsh)-1) |> acf() • Shuffled plot shows no serial correlation

• These two plots in combination show strong effect of ordering

### Measuring effect/importance of ordering for volatility clustering

tr <- wilsh |> logret_1()
shuffled_tr <- tr
shuffled_s <- tr |> as.vector() |> sample(size=nrow(tr))
shuffled_tr$TR <- shuffled_s tr |> get_df() |> ggplot(mapping = aes(x = date, y=TR)) + geom_line() shuffled_tr |> get_df() |> ggplot(mapping = aes(x = date, y=TR)) + geom_line() tr |> abs() |> acf() shuffled_tr |> abs() |> acf() Actual TR Permutation TR Actual TR ACF Permutation TR ACF ## GARCH model ### GARCH (1,1) - normal model How to use rugarch package? • Specify the GARCH model using the ‘ugrachspec’ function • Then specify • variance.model • mean.model • distribution.model logret <- wilsh |> logret_1() garch.N <- ugarchspec(variance.model=list(model='sGARCH', garchOrder=c(1,1)), mean.model = list(armaOrder = c(0,0), include.mean=TRUE), distribution.model = "norm"); garch.N  *---------------------------------* * GARCH Model Spec * *---------------------------------* Conditional Variance Dynamics ------------------------------------ GARCH Model : sGARCH(1,1) Variance Targeting : FALSE Conditional Mean Dynamics ------------------------------------ Mean Model : ARFIMA(0,0,0) Include Mean : TRUE GARCH-in-Mean : FALSE Conditional Distribution ------------------------------------ Distribution : norm Includes Skew : FALSE Includes Shape : FALSE Includes Lambda : FALSE  fit.garch.N <- ugarchfit(spec = garch.N, data=logret); fit.garch.N  *---------------------------------* * GARCH Model Fit * *---------------------------------* Conditional Variance Dynamics ----------------------------------- GARCH Model : sGARCH(1,1) Mean Model : ARFIMA(0,0,0) Distribution : norm Optimal Parameters ------------------------------------ Estimate Std. Error t value Pr(>|t|) mu 0.000678 0.000079 8.5856 0.000000 omega 0.000002 0.000001 2.8917 0.003832 alpha1 0.090834 0.007943 11.4364 0.000000 beta1 0.893290 0.008580 104.1170 0.000000 Robust Standard Errors: Estimate Std. Error t value Pr(>|t|) mu 0.000678 0.000086 7.89699 0.00000 omega 0.000002 0.000006 0.30107 0.76336 alpha1 0.090834 0.060970 1.48981 0.13627 beta1 0.893290 0.071761 12.44813 0.00000 LogLikelihood : 31772.43 Information Criteria ------------------------------------ Akaike -6.6302 Bayes -6.6272 Shibata -6.6302 Hannan-Quinn -6.6291 Weighted Ljung-Box Test on Standardized Residuals ------------------------------------ statistic p-value Lag 24.81 6.342e-07 Lag[2*(p+q)+(p+q)-1] 24.81 2.438e-07 Lag[4*(p+q)+(p+q)-1] 27.70 1.402e-07 d.o.f=0 H0 : No serial correlation Weighted Ljung-Box Test on Standardized Squared Residuals ------------------------------------ statistic p-value Lag 0.3101 0.5776 Lag[2*(p+q)+(p+q)-1] 2.8686 0.4316 Lag[4*(p+q)+(p+q)-1] 3.7794 0.6264 d.o.f=2 Weighted ARCH LM Tests ------------------------------------ Statistic Shape Scale P-Value ARCH Lag 0.01043 0.500 2.000 0.9186 ARCH Lag 0.01799 1.440 1.667 0.9989 ARCH Lag 0.75403 2.315 1.543 0.9500 Nyblom stability test ------------------------------------ Joint Statistic: 179.7119 Individual Statistics: mu 0.06819 omega 22.32490 alpha1 0.13086 beta1 0.19771 Asymptotic Critical Values (10% 5% 1%) Joint Statistic: 1.07 1.24 1.6 Individual Statistic: 0.35 0.47 0.75 Sign Bias Test ------------------------------------ t-value prob sig Sign Bias 2.182 2.915e-02 ** Negative Sign Bias 2.016 4.384e-02 ** Positive Sign Bias 2.799 5.130e-03 *** Joint Effect 48.248 1.886e-10 *** Adjusted Pearson Goodness-of-Fit Test: ------------------------------------ group statistic p-value(g-1) 1 20 478.2 2.065e-89 2 30 418.6 1.241e-70 3 40 813.1 1.962e-145 4 50 782.4 8.811e-133 Elapsed time : 0.4649279  save1 <- cbind(logret, fit.garch.N@fit$sigma, fit.garch.N@fit$z) names(save1) <- c('logret','s','z') save1  logret s z 1980-01-02 -0.0212773984 0.010722744 -2.04757252 1980-01-04 0.0106952891 0.012172801 0.82290777 1980-01-07 0.0053050522 0.011965039 0.38669792 1980-01-08 0.0209431738 0.011467887 1.76710648 1980-01-09 0.0000000000 0.012508564 -0.05421857 1980-01-10 0.0103093697 0.011895049 0.80967916 1980-01-11 0.0000000000 0.011683400 -0.05804787 1980-01-14 0.0051151007 0.011120254 0.39899305 1980-01-15 0.0050890695 0.010674035 0.41323391 1980-01-16 0.0000000000 0.010258012 -0.06611383 ... 2017-12-15 0.0091079755 0.005269317 1.59978602 2017-12-18 0.0062944043 0.005739355 0.97854343 2017-12-19 -0.0037072899 0.005828619 -0.75240576 2017-12-20 -0.0004037957 0.005811807 -0.18617137 2017-12-21 0.0024203320 0.005653465 0.30815359 2017-12-22 -0.0005642204 0.005523522 -0.22493200 2017-12-26 -0.0007258943 0.005392257 -0.26039017 2017-12-27 0.0008065167 0.005275921 0.02432186 2017-12-28 0.0023352267 0.005152582 0.32159219 2017-12-29 -0.0053225932 0.005064389 -1.18489903 glue("{fit.garch.N@model$modeldesc$vmodel}-{fit.garch.N@model$modeldesc$distno}-{fit.garch.N@model$modeldesc$distribution}") sGARCH-1-norm g <- save1$z
g |> as.vector() |> mean() |> round(digits=4)
 -0.0309
g |> as.vector() |> sd() |> round(digits=4)
 1.0006
g |> as.vector() |> skewness() |> round(digits=2)
 -0.59
g |> as.vector() |> kurtosis() |> round(digits=2)
 6.44
g |> as.vector() |> jarque.test()

Jarque-Bera Normality Test

data:  as.vector(g)
JB = 5292.3, p-value < 2.2e-16
alternative hypothesis: greater
g |> acf()
g |> abs() |> acf()
g <- save1$logret g |> as.vector() |> mean() |> round(digits=4)  4e-04 g |> as.vector() |> sd() |> round(digits=4)  0.0107 g |> as.vector() |> skewness() |> round(digits=2)  -0.91 g |> as.vector() |> kurtosis() |> round(digits=2)  21.8 g |> as.vector() |> jarque.test()  Jarque-Bera Normality Test data: as.vector(g) JB = 142514, p-value < 2.2e-16 alternative hypothesis: greater t <- g |> as.vector() |> jarque.test() t$statistic
      JB
142514.5 
t$p.value  0 t$alternative
 "greater"
Interpretation of Garch(1,1) - normal model
• For a normal distribution to be true fitted e(t) should have

• mean=0, sd=1, skewness=0, kurtosis=3

• Jacque bera test should be able to confirm normality

• Kurtosis of our data is around 21.8 . Our fitted distribution is around 6.44 and normal distribution expects it to be 3 => GARCH model is able to explain quite a lot but not all of the heavy tail

• ACF curves for fitted value indicated no serial correlation

• ACF curves for absolute fitted value indicate no serial correlation and thus, no volatility clustering => Variance equation in GARCH model is able to explain all volatility clustering in our data

• Conclusion: GARCH(1,1) - normal model can explain

• Volatility clustering of our data

• Some but not all of heavy tails of the data

• To have better explanation of distribution we need to change distribution equation.

get_fit_params_garch  <- function(s, spec=NULL, solver="solnp"){
if(is.null(spec)){
spec <- ugarchspec(variance.model=list(model='sGARCH', garchOrder=c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean=TRUE),
distribution.model = "std")
}
fit.garch.dist <- ugarchfit(spec = spec, data=s, solver=solver)
return(fit.garch.dist)
}
s <- wilsh |> logret_1(); s |> head()
                     TR
1980-01-02 -0.021277398
1980-01-04  0.010695289
1980-01-07  0.005305052
1980-01-08  0.020943174
1980-01-09  0.000000000
1980-01-10  0.010309370
fit.garch <- s |> get_fit_params_garch(); s |> head()
                     TR
1980-01-02 -0.021277398
1980-01-04  0.010695289
1980-01-07  0.005305052
1980-01-08  0.020943174
1980-01-09  0.000000000
1980-01-10  0.010309370

### Explainability_dist_model

explainability_dist_model <- function(s, spec=NULL, solver="solnp"){
fit.garch <- s |> get_fit_params_garch(spec=spec, solver=solver)
save <- cbind(s, fit.garch@fit$sigma, fit.garch@fit$z)
names(save) <- c('logret','s', 'z')
models <- c("Norm", "Actual",
glue("{fit.garch@model$modeldesc$vmodel}-{fit.garch@model$modeldesc$distno}-{fit.garch@model$modeldesc$distribution}"))

means <- c(0, save$logret |> mean(), save$z |> mean())
sds <- c(1, save$logret |> sd(), save$z |> sd())
skewnesss  <- c(0, save$logret |> skewness(), save$z |> skewness())
kurtosiss <- c(3, save$logret |> kurtosis(), save$z |> kurtosis())
jb_actual  <- save$logret |> as.vector() |> jarque.test() jb_dist <- save$z |> as.vector() |> jarque.test()
JBSTAT.stat <- c("", jb_actual$statistic, jb_dist$statistic )
JBSTAT.pvalue <- c("", jb_actual$p.value, jb_dist$p.value )
JBSTAT.alternative <- c("", jb_actual$alternative, jb_dist$alternative)
df <- data.frame( mean = means,
sd = sds,
skewness = skewnesss,
kurtosis = kurtosiss,
JBSTAT.stat = JBSTAT.stat,
JBSTAT.pvalue = JBSTAT.pvalue,
JBSTAT.alternative = JBSTAT.alternative)
row.names(df) <- models

out <- list(save = save, df=df, parm=fit.garch@fit$coef) return(out) } out <- wilsh |> logret_1() |> explainability_dist_model() out$parm
          mu        omega       alpha1        beta1        shape
7.906871e-04 7.955355e-07 7.125105e-02 9.228052e-01 6.140006e+00 
out$df |> format(digits=4) fitted <- out$save$z fitted |> acf() fitted |> abs() |> acf() actual<- out$save$logret actual |> acf() actual |> abs() |> acf() ACF of fitted e(t) ACF of |fitted e(t)| ACF of actual e(t)-Series Correlation ACF of |actual e(t)|-Volatility Clustering df <- data.frame(date=index(out$save),
s=as.numeric(out$save$s),
TR=as.numeric(out$save$logret))
ggplot(data=df) + geom_line(mapping = aes(x=date, y=s))
ggplot(data=df) + geom_line(mapping = aes(x=date, y=TR))

## Exercise 12

load("W1_Exercise2_FRED_gold.gz")
g <- gold |> filterSeries() |> logret_1()
g |> head(2)
                    TR
1980-01-03  0.12500543
1980-01-04 -0.07532201
g |> tail(2)
                    TR
2017-12-27 0.011674890
2017-12-28 0.009025894
out <- g |> explainability_dist_model()
out$parm |> round(digits=6)  mu omega alpha1 beta1 shape 0.000020 0.000000 0.058965 0.940035 4.705043  out$df |> format(digits=4)
fitted <-  out$save$z
fitted |> acf()
fitted |> abs() |> acf()
actual<-  out$save$logret
actual |> acf()
actual |> abs() |> acf()

## VaR and ES for GARCH Bootstrap

get_garch_series_from_s <- function(s, n=100000, fitpars=NULL, set_seed=TRUE, solver="solnp"){

if (set_seed){
RNGkind(sample.kind="Rounding")
set.seed(123789)
}
if(is.null(fitpars)){
fitpars <- s |> get_fit_params_garch(solver = solver)
}

boot.garch  <- ugarchboot(fitpars,
method = "Partial",
sampling="raw",
n.bootpred = n,
solver="solnp"
)
return(boot.garch@fseries |> as.vector())
}

### VaR and ES nextday (Jan2, 2018) with 95% CI

logret <- wilsh |> logret_1()
s  <- logret |> get_garch_series_from_s()

wilsh |> head(1) |> index()
 "1979-12-31"
wilsh |> tail(1) |> index()
 "2017-12-29"
VaR <-  s |> get_sim_VaR(alpha=0.05); VaR
 -0.007456152
ES <- s |> get_sim_ES(VaR=VaR);ES
 -0.01094603

### VaR and ES nextday (Oct19, 1987) with 95% CI

logret <- wilsh |> logret_1() |> filterSeries(date_filter = '1979-12-31/1987-10-19')
s  <- logret |> get_garch_series_from_s()

logret |> head(1) |> index()
 "1980-01-02"
logret |> tail(1) |> index()
 "1987-10-19"
VaR <-  s |> get_sim_VaR(alpha=0.05); VaR
 -0.07382125
ES <- s |> get_sim_ES(VaR=VaR);ES
 -0.1058245

### GARCH 1day many input model

• 2008-09-15 -Lehmann failure

• 1987-10-19 - Largest single day decline in US stock market

do_analysis_garch_nextday <- function(rvec,
date.start="1979-12-31",
date.ends=c("2017-12-29"),
fitspar=NULL,
n=100000,
alpha=0.05,
solver="solnp"){

VaRs <- list()
ESs <- list()
for(i in date.ends){
filter_str = glue('{date.start}/{i}')
rvec_filter <- rvec |> filterSeries(date_filter=filter_str)
s <- rvec_filter |> get_garch_series_from_s(fitpars=fitspar,
n=n,
solver=solver)
VaR <-  s |> get_sim_VaR(alpha=alpha)
ES <- s |> get_sim_ES(VaR=VaR)
VaRs <- c(VaRs, VaR)
ESs <- c(ESs,ES)
}

return(data.frame(day=date.ends |> as.Date(),
Var=VaRs |> as.double(),
ES=ESs |> as.double()))
}

wilsh |> logret_1() |> do_analysis_garch_nextday(date.ends=c("2017-12-29","2008-09-15","1987-10-19")) |> format(digits=4)

### Rolling 1-day VaR

s <- wilsh |> logret_1()
spec <- ugarchspec(variance.model=list(model='sGARCH', garchOrder=c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean=TRUE),
distribution.model = "std")

s2016 <- s |> filterSeries(date_filter = '1980-01-01/2016-12-31')
n2016 <- length(s2016); n2016
 9332
roll.garch <- ugarchroll(spec=spec,
data=s,
start=1,
n.start =n2016,
forecast.length=1,
refit.every = 100,
refit.window = "recursive",
calculate.VaR = TRUE,
VaR.alpha = 0.05,
keep.coef = TRUE)
roll.garch@forecast$VaR |> head(1) roll.garch@forecast$VaR |> tail(1)
s_from2017 <- s |> filterSeries(date_filter = '2016-12-31/2018-01-01')
s_from2017 |> head(1)
                    TR
2017-01-03 0.008281823
s_from2017 |> tail(1)
                     TR
2017-12-29 -0.005322593
roll_df <- data.frame(date=index(s_from2017), TR=s_from2017$TR, VaR=roll.garch@forecast$VaR)

colnames(roll_df) <- c('date', 'TR', 'VaR')
roll_df
sum(roll_df$TR < roll_df$VaR)*100/nrow(roll_df) |> round(digits=1)
 2.390438
roll.garch@forecast$VaR |> head(5) |> as.vector()  $alpha(5%)
 -0.008529988 -0.008837648 -0.009013061 -0.008797251 -0.008579853
ggplot(data=roll_df)+
geom_line(mapping = aes(x=date, y=TR, color='TR')) +
geom_line(mapping = aes(x=date, y=VaR, color='VaR')) get_garch_VaR_plot <- function(s,
spec=NULL,
alpha=0.05,
refit.every=1,
date.start='1979-12-31',
date.fit='2016-12-31',
date.predict='2018-01-01'){
if(is.null(spec)){
spec = ugarchspec(variance.model=list(model='sGARCH', garchOrder=c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean=TRUE),
distribution.model = "std")
}

sfit <- s |> filterSeries(date_filter = glue('{date.start}/{date.fit}'))
nfit <- length(sfit)
roll.garch <- ugarchroll(spec=spec,
data=s,
start=1,
n.start=nfit,
forecast.length=1,
refit.every = refit.every,
refit.window = "recursive",
calculate.VaR = TRUE,
VaR.alpha = alpha,
keep.coef = TRUE)

spredict <- s |> filterSeries(date_filter = glue('{date.fit}/{date.predict}'))
roll_df <- data.frame(date=index(spredict),
TR=spredict$TR, VaR=roll.garch@forecast$VaR)
colnames(roll_df) <- c('date', 'TR', 'VaR')

p <- ggplot(data=roll_df)+
geom_line(mapping = aes(x=date, y=TR, color='TR')) +
geom_line(mapping = aes(x=date, y=VaR, color='VaR')) +
labs(title=glue("1 day {(1-alpha)*100}% VaR from {date.fit} to {date.predict}" ), subtitle = glue("TR < VaR for {sum(roll_df$TR < roll_df$VaR)*100/nrow(roll_df)}%") )
return(p)

}
wilsh |> logret_1() |> get_garch_VaR_plot(refit.every = 100, alpha = 0.05) ## Exercise 13

g <- gold |> filterSeries() |> logret_1()
g |> head(2)
                    TR
1980-01-03  0.12500543
1980-01-04 -0.07532201
g |> tail(2)
                    TR
2017-12-27 0.011674890
2017-12-28 0.009025894
g |> do_analysis_garch_nextday(date.ends=c("2017-12-29","2008-09-15","1987-10-19")) |> format(digits=5)
Warning in RNGkind(sample.kind = "Rounding"): non-uniform 'Rounding' sampler
used

Warning in RNGkind(sample.kind = "Rounding"): non-uniform 'Rounding' sampler
used

Warning in RNGkind(sample.kind = "Rounding"): non-uniform 'Rounding' sampler
used
g |> get_garch_VaR_plot(refit.every = 100, alpha = 0.05) ## Quiz 1

us2yen <-getSeries("DEXJPUS")
yen2us <- 1/us2yen
yen2us |> head(2)
                    TR
1979-12-31 0.004161465
1980-01-02 0.004193751
yen2us |> tail(2)
                    TR
2017-12-28 0.008861320
2017-12-29 0.008873902
tr <- yen2us |> logret_1()
tr |> acf() tr <- yen2us |> logret_1()
shuffled_tr <- tr
shuffled_s <- tr |> as.vector() |> sample(size=nrow(tr))
shuffled_tr$TR <- shuffled_s tr |> get_df() |> ggplot(mapping = aes(x = date, y=TR)) + geom_line() shuffled_tr |> get_df() |> ggplot(mapping = aes(x = date, y=TR)) + geom_line() tr |> abs() |> acf() shuffled_tr |> abs() |> acf() s <- yen2us |> logret_1() spec <- ugarchspec(variance.model=list(model='sGARCH', garchOrder=c(1,1)), mean.model = list(armaOrder = c(0,0), include.mean=TRUE), distribution.model = "std") fit.garch <- ugarchfit(spec = spec, data=s, solver = 'hybrid', solver.control = list(), fit.control = list(stationarity = 1, fixed.se = 0, scale = 0, rec.init = 'all', trunclag = 1000) ) save <- cbind(s, fit.garch@fit$sigma, fit.garch@fit$z) names(save) <- c('logret','s', 'z') fit.garch@fit$coef
           mu         omega        alpha1         beta1         shape
-9.894718e-05  1.079617e-07  3.294518e-02  9.661616e-01  4.960222e+00 
out <- yen2us |> logret_1() |> explainability_dist_model(solver="hybrid")
out$parm |> round(digits=6)  mu omega alpha1 beta1 shape -0.000099 0.000000 0.032945 0.966162 4.960222  yen2us |> logret_1() |> do_analysis_garch_nextday(date.ends=c("2018-01-01","2008-09-15","1987-10-19"), solver='hybrid') |> format(digits=4) Warning in RNGkind(sample.kind = "Rounding"): non-uniform 'Rounding' sampler used Warning in RNGkind(sample.kind = "Rounding"): non-uniform 'Rounding' sampler used Warning in RNGkind(sample.kind = "Rounding"): non-uniform 'Rounding' sampler used yen2us |> logret_1() |> tail(1)  TR 2017-12-29 0.001418817 yen2us |> logret_1() |> do_analysis_garch_nextday(date.ends=c("2017-12-29"), solver='hybrid') |> format(digits=4) Warning in RNGkind(sample.kind = "Rounding"): non-uniform 'Rounding' sampler used ## Quiz 2 bp2us <-getSeries("DEXUSUK") bp2us |> head(2)  TR 1979-12-31 2.2185 1980-01-02 2.2427 bp2us |> tail(2)  TR 2017-12-28 1.3450 2017-12-29 1.3529 # yen2us <- 1/us2yen # yen2us |> head(2) # yen2us |> tail(2) s <- bp2us |> logret_1() s |> head(2)  TR 1980-01-02 0.010849205 1980-01-03 -0.001204631 s |> tail(2)  TR 2017-12-28 0.003500544 2017-12-29 0.005856424 tr <- bp2us |> logret_1() tr |> acf() tr <- bp2us |> logret_1() shuffled_tr <- tr shuffled_s <- tr |> as.vector() |> sample(size=nrow(tr)) shuffled_tr$TR <- shuffled_s

tr |>
get_df() |>
ggplot(mapping = aes(x = date, y=TR)) +
geom_line() shuffled_tr |>
get_df() |>
ggplot(mapping = aes(x = date, y=TR)) +
geom_line() tr |> abs() |> acf() shuffled_tr |> abs() |> acf() out <- bp2us |> logret_1() |> explainability_dist_model()
out\$parm |> round(digits=6)
      mu    omega   alpha1    beta1    shape
0.000063 0.000000 0.044625 0.949961 6.980541 
bp2us |> logret_1() |> do_analysis_garch_nextday(date.ends=c("2017-12-29","2008-09-15","1987-10-19")) |> format(digits=5)