Package {stochvolTMB}

Backtest Value at Risk and Expected Shortfall

Description

Tests if the exceedances of the Value at Risk (returns below -var) have the expected frequency and are independent over time.

Usage

backtest_var(x)

Arguments

Details

Let I_t = 1 if y_t < -VaR_t and p = 1 - \alpha the expected exceedance rate.

• Unconditional coverage (Kupiec 1995): a likelihood ratio test of P(I_t = 1) = p. A small p-value means too many or too few exceedances.

• Independence (Christoffersen 1998): a likelihood ratio test that I_t is a Markov chain with P(I_t = 1 \mid I_{t-1} = 0) = P(I_t = 1 \mid I_{t-1} = 1). A small p-value means that exceedances come in clusters, which happens when the model reacts too slowly to changes in volatility.

• Conditional coverage: the sum of the two statistics, with two degrees of freedom.

es_ratio is the mean of -y_t / ES_t over the exceedances. It is close to one if the Expected Shortfall is correct, and larger than one if the model underestimates the size of the losses in the tail.

Value

A data.table with one row for each level and the columns level, n, expected and observed (number of exceedances), rate (observed rate), lr_uc, p_uc, lr_ind, p_ind, lr_cc, p_cc and es_ratio.

References

Kupiec, P. H. (1995). Techniques for verifying the accuracy of risk measurement models. The Journal of Derivatives, 3(2), 73-84. doi:10.3905/jod.1995.407942

Christoffersen, P. F. (1998). Evaluating interval forecasts. International Economic Review, 39(4), 841-862. doi:10.2307/2527341

See Also

var_es

Extract parameter estimates and their covariance matrix

Description

Extract parameter estimates and their covariance matrix

Usage

## S3 method for class 'stochvolTMB'
coef(object, scale = c("natural", "estimation"), ...)

## S3 method for class 'stochvolTMB'
vcov(object, scale = c("natural", "estimation"), ...)

Arguments

Value

coef returns a named vector. vcov returns a named matrix.

Examples

opt <- estimate_parameters(spy$log_return, model = "leverage", silent = TRUE)
coef(opt)
sqrt(diag(vcov(opt)))
coef(opt, scale = "estimation")

Fit several models and compare them

Description

Fits each model with estimate_parameters and returns a table ranked by AIC.

Usage

compare_models(
  data,
  models = c("gaussian", "t", "leverage", "skew_gaussian", "skew_gaussian_leverage",
    "skew_t", "skew_t_leverage"),
  ...
)

Arguments

Details

The table also shows if the optimizer converged and if the Hessian is positive definite. A model that fails one of these checks has unreliable estimates and standard errors, even when its AIC is small. A model that gives an error is reported with NA values and the error message.

Value

A data.table with one row for each model, ordered by AIC, and the columns model, n_par, log_lik, aic, bic, delta_aic (the difference to the smallest AIC), converged, pd_hessian and message. The fitted models are in the attribute "fits", a named list of stochvolTMB objects.

Examples

comparison <- compare_models(spy$log_return, models = c("gaussian", "t", "leverage"))
comparison
best <- attr(comparison, "fits")[[comparison$model[1]]]

Run the shiny demo (deprecated)

Description

Use run_shiny_demo instead. Without arguments, demo() runs the shiny demo with a deprecation warning. With arguments, it calls demo, so demo(topic) still works when stochvolTMB is attached.

Usage

demo(...)

Arguments

Value

Called for its side effect.

Estimate parameters for the stochastic volatility model

Description

Estimate parameters of a stochastic volatility model with a latent log-volatility following an autoregressive process of order one with normally distributed noise. The following distributions are implemented for the observed process:

• "gaussian": Gaussian distribution.

• "t": t-distribution, scaled to unit variance. Estimates the degrees of freedom df.

• "leverage": Gaussian distribution where the noise of the latent process is correlated with the observational noise. Estimates the correlation rho.

• "skew_gaussian": Skew gaussian distribution with mean zero and unit variance. Estimates the skewness of the observational noise.

• "skew_gaussian_leverage": Skew gaussian distribution with leverage. The skew gaussian return shock z_t drives the volatility shock, eta_t = rho * z_t + sqrt(1 - rho^2) * e_t, where e_t is standard gaussian. Estimates both the skewness and rho.

• "skew_t": Skew t distribution of Jones and Faddy (2003), scaled to mean zero and unit variance. The left and the right tail have their own degrees of freedom, df_left and df_right, so the model can have heavier tails for losses than for gains. With equal values it is the "t" model.

• "skew_t_leverage": Skew t distribution with leverage. As in "skew_gaussian_leverage", the return shock z_t drives the volatility shock, eta_t = rho * z_t + sqrt(1 - rho^2) * e_t. Estimates df_left, df_right and rho.

The parameters is estimated by minimizing the negative log-likelihood (nll) and the latent log-volatility is integrated out by applying the Laplace approximation.

Usage

estimate_parameters(data, model = "gaussian", opt.control = NULL, ...)

Arguments

Value

Object of class stochvolTMB

Skew gaussian models

The skew gaussian models report two numbers for the asymmetry of the standardized noise z:

• alpha: the shape parameter of the skew normal density. alpha = 0 gives the gaussian distribution and alpha -> Inf the half-gaussian distribution. It has no direct statistical meaning.

• skewness: the skewness E[z^3] of the noise, a monotone function of alpha. With delta = alpha / sqrt(1 + alpha^2) and mu = delta * sqrt(2 / pi), skewness = (4 - pi) / 2 * mu^3 / (1 - mu^2)^(3 / 2). For example, alpha = 0.5, 1, 2 and 5 give skewness = 0.02, 0.14, 0.45 and 0.85.

The models estimate the skewness, not alpha. Close to zero the skewness grows like alpha^3, so the likelihood is flat in alpha at alpha = 0 and the optimizer does not move from there.

Use summary(object, report = "transformed") to get skewness and alpha with standard errors. Use skewness for inference. Use alpha when you need the shape parameter, for example in sim_sv or sn::dsn. Close to zero skewness, alpha changes fast with the skewness, so the standard error of alpha is very large even when the skewness is well determined. In the "skew_gaussian_leverage" model the skewness is that of the return shock, not of the returns.

Skew t models

The left tail of the noise decays like a Student t distribution with df_left degrees of freedom and the right tail like one with df_right degrees of freedom. A smaller value means a heavier tail. If df_left < df_right, large losses are more likely than large gains of the same size, and the noise is skewed to the left. Both values are greater than 2, so that the variance exists. They are estimated on the scale log(df - 2).

A light tail is hard to tell apart from a gaussian tail. If one tail is light, its estimated degrees of freedom can be very large and have a large standard error, and nlminb can report false convergence. This means that the data show no heavy tail on that side. The same is true for df in the "t" model. In both models the degrees of freedom are capped at 10,000, because the density cannot be calculated accurately for larger values. Short series (a few hundred observations) can also confuse heavy tails with a very noisy log-volatility.

Hansen's (1994) skewed t distribution is not used, because its log density has a jump in the second derivative at the mode. This makes the Laplace approximation inaccurate and the optimization unstable.

References

Jones, M. C. and Faddy, M. J. (2003). A skew extension of the t-distribution, with applications. Journal of the Royal Statistical Society: Series B, 65(1), 159-174. doi:10.1111/1467-9868.00378

Examples

# load data
data("spy")

# estimate parameters 
opt <- estimate_parameters(spy$log_return, model = "gaussian")

# get parameter estimates with standard error
estimates <- summary(opt)

# plot estimated volatility with 95 % confidence interval
plot(opt, include_ci = TRUE)

# skew gaussian model with leverage
opt_sgl <- estimate_parameters(spy$log_return, model = "skew_gaussian_leverage", silent = TRUE)
summary(opt_sgl, report = "transformed")

Value at Risk and Expected Shortfall for several days ahead

Description

Simulates future returns with predict.stochvolTMB and calculates the Value at Risk (VaR) and Expected Shortfall (ES) for each day ahead, both of the return on that day and of the cumulative return (the sum of the log-returns) up to that day. The cumulative values for step = 10 are, for example, the 10-day VaR and ES.

Usage

forecast_var_es(
  object,
  steps = 10L,
  level = c(0.95, 0.99),
  nsim = 100000L,
  include_parameters = TRUE
)

Arguments

Details

VaR and ES are reported as positive numbers when they are losses: if q is the 1 - \alpha quantile of the simulated returns, then VaR = -q and ES = -E[y \mid y \le q]. By default, the simulation includes the uncertainty of the parameters. The Monte Carlo error decreases with nsim. For one day ahead, var_es is exact and faster.

Value

A data.table with the columns step, level, var, es, var_cumulative and es_cumulative.

See Also

var_es

Examples

opt <- estimate_parameters(spy$log_return, model = "skew_t_leverage", silent = TRUE)
set.seed(1)
forecast_var_es(opt, steps = 10, level = 0.99)

Construct objective function with derivatives using MakeADFun

Description

Construct objective function with derivatives using MakeADFun

Usage

get_nll(data, model = "gaussian", ...)

Arguments

Value

List of components returned from MakeADFun.

Filter on a grid for the log-volatility

Description

Keeps the predictive density of h_t on a grid. For each observation it calculates the predictive distribution function of y_t, updates the density with y_t and moves it one step forward. Grid points with negligible weight are skipped, so the cost per step is small.

Usage

grid_filter(object, grid_size = 500, width = 8)

Arguments

Value

A list with prob, the predictive distribution function at each y_t; weights, a list of length n + 1 where element t holds the grid indices idx and the normalized predictive weights w of h_t given y_1, ..., y_{t-1} (element n + 1 is the forecast); and scale, the conditional standard deviation sigma_y * exp(h / 2) at each grid point.

Logit transformation from the real line to (-1, 1).

Description

Logit transformation from the real line to (-1, 1).

Usage

logit(x)

Arguments

Value

double

Number of loaded DLLs that contain a TMB model

Description

Number of loaded DLLs that contain a TMB model

Usage

n_tmb_dlls()

Value

Integer.

Exact one-step-ahead residuals by filtering on a grid for the log-volatility

Description

Exact one-step-ahead residuals by filtering on a grid for the log-volatility

Usage

osa_grid(object, grid_size = 500)

Arguments

Value

Vector of one-step-ahead residuals for all observations.

Plot the estimated latent volatility process

Description

Displays the estimated latent volatility process over time.

Usage

## S3 method for class 'stochvolTMB'
plot(x, ..., include_ci = TRUE, plot_log = TRUE, dates = NULL, forecast = NULL)

Arguments

Value

ggplot object with plot of estimated estimated volatility.

Diagnostic plots based on one-step-ahead residuals

Description

Shows three panels that check if the fitted model describes the data:

• A Q-Q plot of the one-step-ahead (OSA) residuals against the standard gaussian distribution. If the model is correct, the points follow the diagonal. Bent tails show that the error distribution is wrong, for example that the returns have heavier tails or more skewness than the model allows.

• The autocorrelation function (ACF) of the residuals. Bars outside the dashed lines show dependence in the returns that the model does not capture.

• The ACF of the squared residuals. Bars outside the dashed lines show volatility clustering that the model does not capture.

Usage

plot_diagnostics(x, res = NULL, lag_max = 20, ...)

Arguments

Details

By default, the residuals are calculated for all observations with the exact grid method of residuals, which takes a few seconds. To reuse residuals, calculate them once with residuals and give them in res.

Value

ggplot object with the three diagnostic panels.

Examples

data(spy)
opt <- estimate_parameters(spy$log_return, model = "t", silent = TRUE)
plot_diagnostics(opt)

Add predicted volatility.

Description

Adds predicted volatility to the volatility plot.

Usage

plot_forecast(p, forecast, include_ci = TRUE)

Arguments

Plot returns with bands from the fitted model

Description

Plots the observed returns together with the central level interval of the fitted conditional distribution of each return, given the estimated log-volatility h_t. The interval uses the error distribution of the model, so it is asymmetric for the skew gaussian models and wider for the t model. Returns outside the interval are marked in red.

Usage

plot_returns(x, level = 0.95, dates = NULL)

Arguments

Details

The subtitle shows the share of returns outside the interval. For a good model it is close to 1 - level. The bands use the smoothed estimate of h_t, which also uses y_t. The share is therefore an optimistic, in-sample check. Use plot_diagnostics for a check based on one-step-ahead residuals.

Value

ggplot object.

Examples

data(spy)
opt <- estimate_parameters(spy$log_return, model = "skew_gaussian_leverage", silent = TRUE)
plot_returns(opt, dates = spy$date)

Plot returns with the one-step-ahead Value at Risk

Description

Plots the returns with -VaR and -ES from var_es and marks the exceedances, the returns below -VaR.

Usage

plot_var(x, level = NULL, dates = NULL)

Arguments

Value

ggplot object.

See Also

var_es, backtest_var

Examples

opt <- estimate_parameters(spy$log_return, model = "skew_t_leverage", silent = TRUE)
plot_var(var_es(opt, level = 0.99), dates = spy$date)

Predict future returns and future volatilities

Description

Takes a stochvolTMB object and produces draws from the predictive distribution of the latent volatility and future log-returns.

Usage

## S3 method for class 'stochvolTMB'
predict(object, steps = 1L, nsim = 10000, include_parameters = TRUE, ...)

Arguments

Value

List of simulated values from the predictive distribution of the latent volatilities and log-returns.

Calculate one-step-ahead (OSA) residuals for stochastic volatility model.

Description

The OSA residual of observation y_t is qnorm(F_t(y_t)), where F_t is the predictive distribution function of y_t given y_1, ..., y_{t-1} at the estimated parameters. If the model is correct, the residuals are independent and standard gaussian.

Usage

## S3 method for class 'stochvolTMB'
residuals(
  object,
  conditional = NULL,
  method = c("grid", "tmb"),
  grid_size = 500,
  ...
)

Arguments

Details

Two methods are available:

• "grid" (default): an exact filter on a grid for the log-volatility h_t. Because h_t is one-dimensional, the predictive distribution can be calculated by summation over the grid, without the Laplace approximation. It calculates the residuals for all observations in a few seconds.

• "tmb": oneStepPredict with method = "oneStepGeneric". It uses the Laplace approximation and takes about one to three seconds per observation, so by default only the last 100 observations are used. It is not available for the "skew_t_leverage" model, where y_t also enters the model for h_{t+1}.

Value

Vector of one-step-ahead residuals. If the model is correctly specified, these should be standard normal.

References

Thygesen, U. H., Albertsen, C. M., Berg, C. W., Kristensen, K. and Nielsen, A. (2017). Validation of ecological state space models using the Laplace approximation. Environmental and Ecological Statistics, 24(2), 317-339. doi:10.1007/s10651-017-0372-4

Draw from the skew gaussian distribution standardized to mean zero and unit variance

Description

Draw from the skew gaussian distribution standardized to mean zero and unit variance

Usage

rsn_standard(n, alpha)

Arguments

Value

Numeric vector of length n.

Run the shiny demo

Description

Starts a shiny application where you can simulate from the models and see the effect of the parameters. The shiny package must be installed.

Usage

run_shiny_demo()

Value

Called for its side effect.

Simulate log-returns from a stochastic volatility model

Description

This function draws the initial log-volatility (h_t) from its stationary distribution, meaning that h_0 is drawn from a gaussian distribution with mean zero and standard deviation sigma_h / sqrt(1 - phi^2). h_{t+1} is then simulated from its conditional distribution given h_t, which is N(phi*h_t, sigma_h). Log-returns (y_t) is simulated from its conditional distribution given the latent process h. If model = "gaussian", then y_t given h_t is gaussian with mean zero and standard deviation equal to sigma_y*exp(h_t / 2). Heavy tail returns can be obtained by simulating from the t-distribution by setting model = "t". How heavy of a tail is specified by the degree of freedom parameter df. Note that the observations are scaled by sqrt((df-2)/2) so that the error term has variance equal to one. Asymmetric returns are obtained from the "skew_gaussian" model. How asymmetric is governed by the skewness parameter alpha. The so called leverage model, where we allow for correlation between log-returns and volatility can be simulated by setting model to "leverage" and specifying the correlation parameter rho. In the "skew_gaussian_leverage" model the return shock z_t is skew gaussian with mean zero, unit variance and shape alpha, and it drives the volatility shock: h_{t+1} = phi * h_t + sigma_h * (rho * z_t + sqrt(1 - rho^2) * e_t), where e_t is standard gaussian.

Usage

sim_sv(
  param = list(phi = 0.9, sigma_y = 0.4, sigma_h = 0.2, df = 4, alpha = -2, rho = -0.7),
  nobs = 1000L,
  seed = NULL,
  model = "gaussian"
)

Arguments

Value

data.table with columns y (observations) and h (latent log-volatility).

Simulate from the asymptotic distribution of the parameter estimates

Description

Sampling is done on the scale the parameters were estimated. The standard deviations are simulated on log-scale and the persistence is simulated on logit scale. The same is true for the correlation parameter in the leverage model.

Usage

simulate_parameters(object, nsim = 1000)

Arguments

Value

matrix of simulated values.

Skew t distribution of Jones and Faddy

Description

Density, distribution function, quantile function and random generation for the skew t distribution of Jones and Faddy (2003), standardized to mean zero and unit variance.

Usage

dskew_t(x, df_left, df_right, log = FALSE)

pskew_t(q, df_left, df_right)

qskew_t(p, df_left, df_right)

rskew_t(n, df_left, df_right)

Arguments

Details

Before standardization, the density with parameters a, b > 0 is

f(t) = C \left(1 + \frac{t}{\sqrt{a + b + t^2}}\right)^{a + 1/2} \left(1 - \frac{t}{\sqrt{a + b + t^2}}\right)^{b + 1/2}, \qquad C = \frac{1}{2^{a + b - 1} B(a, b) \sqrt{a + b}}.

The left tail decays like a Student t with 2a degrees of freedom and the right tail like a Student t with 2b degrees of freedom. The functions use df_left = 2a and df_right = 2b. If they are equal, the distribution is the Student t distribution scaled to unit variance. If df_left < df_right, the left tail is heavier and the distribution is skewed to the left. Both must be greater than 2 for the variance to exist.

The density is smooth everywhere. This is important for the Laplace approximation, which uses the second derivative of the log density.

Value

Numeric vector.

References

Jones, M. C. and Faddy, M. J. (2003). A skew extension of the t-distribution, with applications. Journal of the Royal Statistical Society: Series B, 65(1), 159-174. doi:10.1111/1467-9868.00378

Shape parameter of the skew gaussian distribution with a given skewness

Description

Shape parameter of the skew gaussian distribution with a given skewness

Usage

skewness_to_alpha(skewness)

Arguments

Value

The shape parameter alpha.

Daily closing prices for the S&P500 from 2005 to 2018.

Description

A dataset containing the prices and log-returns of the S&P500 from 2005 to 2018

Usage

spy

Format

A data frame with 3522 rows and 3 variables:

date

date

price

price, in US dollars

log_return

logarithmic return

...

Standardized error distribution of a fitted model

Description

The distribution of y_t / (sigma_y * exp(h_t / 2)) given h_t. It has mean zero and unit variance.

Usage

standard_error_distribution(object)

Arguments

Value

A list of vectorized functions: cdf, logpdf, quantile and partial_mean, where partial_mean(z) is E[Z 1(Z <= z)].

Summary tables of model parameters

Description

Extract parameters, transformed parameters and latent log volatility along with standard error, z-value and p-value

Usage

## S3 method for class 'stochvolTMB'
summary(object, ..., report = c("all", "fixed", "transformed", "random"))

Arguments

Value

data.table with parameter estimates, standard error, z-value and approximated p-value.

Calculate quantiles based on predictions from the predictive distribution

Description

Calculate quantiles based on predictions from the predictive distribution

Usage

## S3 method for class 'stochvolTMB_predict'
summary(object, ..., quantiles = c(0.025, 0.975), predict_mean = TRUE)

Arguments

Value

A list of data.tables. One for y, h and h_exp.

One-step-ahead Value at Risk and Expected Shortfall

Description

Calculates the one-day Value at Risk (VaR) and Expected Shortfall (ES) of the returns from the predictive distribution of y_t given y_1, ..., y_{t-1}, for every observation and for the next day after the data.

Usage

var_es(object, level = c(0.95, 0.99), grid_size = 500)

Arguments

Details

The predictive distribution is calculated exactly with the same filter on a grid for the log-volatility as residuals.stochvolTMB. It is a mixture over h_t of the conditional distribution of the returns. For a level \alpha, for example 0.99, let q_t be the 1 - \alpha quantile of the predictive distribution. Then

VaR_t = -q_t, \qquad ES_t = -E[y_t \mid y_t \le q_t, y_1, \ldots, y_{t-1}].

Both are reported as positive numbers when they are losses, on the scale of the returns. A return is an exceedance if y_t < -VaR_t.

The parameters are fixed at their estimates from the full data, and the uncertainty of the parameters is not included. The values for the observed data are therefore in-sample, and a backtest with backtest_var on them is slightly optimistic. The last row (time = nobs + 1) is the forecast for the next day. For more than one day, use forecast_var_es.

Value

A data.table with one row for each time point and level, and the columns time (1 to nobs + 1), level, y (the observed return, NA for the forecast), var, es and exceedance (y < -var).

See Also

backtest_var, forecast_var_es

Examples

opt <- estimate_parameters(spy$log_return, model = "skew_t_leverage", silent = TRUE)
risk <- var_es(opt, level = c(0.95, 0.99))
# Forecast for the next day
risk[is.na(y)]
backtest_var(risk)