Оценка риска при использовании пакета fGarch в R

Question

Я пытаюсь провести подобный анализ с McNeil & Frey в их статье "Оценка мер риска, связанных с хвостом, для гетероскедастических финансовых временных рядов: подход с экстремальной ценностью", но я столкнулся с проблемой при реализации моделей.

Подход заключается в подборе модели AR(1)-GARCH(1,1) для оценки прогноза VaR на один день с использованием окна из 1000 наблюдений.

Я смоделировал данные, которые должны нормально работать с моей моделью, и я предполагаю, что, если я буду делать это правильно, наблюдаемый уровень покрытия должен быть близок к теоретическому. Однако он всегда ниже теоретического уровня охвата, и я не знаю почему.

Я полагаю, что это, как расчет расчетного VaR делается

VaR_hat = mu_hat + sigma_hat * qnorm(альфа)

, но я могу ошибаться. Я попытался найти похожие вопросы здесь, в стеке, но я не нашел ни одного.

Как я подхожу к этому, можно суммировать в три этапа.

1. Смоделируйте 2000 AR(1)-GARCH(1,1) наблюдений и подберите соответствующую модель и извлеките однодневный прогноз условного среднего и стандартного отклонения, используя окно из 1000 наблюдений (тем самым делая 1000 прогнозов).

2. Используйте прогнозируемые значения и нормальный квантиль для расчета VaR для требуемого уровня достоверности.

3. Проверьте, близок ли уровень покрытия к теоретическому.

Если бы кто-то мог помочь мне, я был бы чрезвычайно благодарен, и если я неясен в моем формировании, пожалуйста, просто скажите мне, и я постараюсь найти лучшее объяснение этой проблемы.

Код, который я использую, прилагается ниже. заранее спасибо

library(fGarch)

nObs <- 2000                  # Number of observations.
quantileLevel <- 0.95         # Since we expect 5% exceedances.  
from <- seq(1,1000)           # Lower index vector for observations in model.
to <- seq(1001,2000)          # Upper index vector for observations in model.
VaR_vec <- rep(0,(nObs-1000)) # Empty vector for storage of 1000 VaR estimates.

# Specs for simulated data (including AR(1) component and all components for GARC(1,1)).
spec = garchSpec(model = list(omega = 1e-6, alpha = 0.08, beta = 0.91, ar = 0.10),
             cond.dist = 'norm')
# Simulate 1000 data points.
data_sim <- c(garchSim(spec, n = nObs, n.start = 1000))

for (i in 1:1000){
  # The rolling window of 1000 observations.
  data_insert <- data_sim[from[i]:to[i]]
  # Fitting an AR(1)-GARCH(1,1) model with normal cond.dist.
  fitted_model <- garchFit(~ arma(1,0) + garch(1,1), data_insert,
                       trace = FALSE,
                       cond.dist = "norm")
  # One day ahead forecast of conditional mean and standard deviation.
  predict(fitted_model, n.ahead = 1)
  prediction_model <- predict(fitted_model, n.ahead = 1)
  mu_pred <- prediction_model$meanForecast
  sigma_pred <- prediction_model$standardDeviation
  # Calculate VaR forecast
  VaR_vec[i] <- mu_pred + sigma_pred*qnorm(quantileLevel)

  if (length(to)-i != 0){
    print(c('Countdown, just',(length(to) - i),'iterations left'))
  } else {
    print(c('Done!'))
  }
}

# Exctract only the estiamtes ralated to the forecasts.
compare_data_sim <- data_sim[1001:length(data_sim)]
hit <- rep(0,length(VaR_vec))
# Count the amount of exceedances. 
for (i in 1:length(VaR_vec)){
  hit[i] <- sum(VaR_vec[i] <= compare_data_sim[i])
}


plot(data_sim[1001:2000], type = 'l', 
 ylab = 'Simulated data', main = 'Illustration of one day ahead prediction of 95%-VaR')
lines(VaR_vec, col = 'red')

cover_prop <- sum(hit)/length(hit)
print(sprintf("Diff theoretical level and VaR coverage = %f", (1-quantileLevel) - cover_prop))