library(tidyverse) knitr::opts_chunk$set(message = FALSE, warning = FALSE, fig.width = 8, fig.height = 4.5, fig.align = 'center', out.width='95%', dpi = 200) # library(ESG) devtools::load_all() theme_set(theme_classic())
Introduction
This package provides functions for generating economic scenario sets for (Re)insurer Economic Capital Models (ECMs).
The motivation is to make seamless functionality for calibration and simulation of financial variables to be used standalone or within (Re)insurer ECMs.
Market variables of interest include:
- Inflation
- Interest Rates (real & nominal)
- Equity Returns
This package relies heavily on the tidyverse.
Core Functions
For ease of reference, functions have been grouped into two broad categories:
- Calibration (Cal*): Given user specified data returns the parameters of a given model. e.g.
CalVasicek1f() - Simulation (Sim*): Given user specified parameters returns simulation output. e.g.
SimVasicek1f()
Datasets
Some datasets are available with the package. Those can be accessed using the data command. Below is a list of the package's datasets:
data(package = "ESG")$results[, "Item"]
Interest Rates
The package includes tools to model various short rate models of the form: $dr_t = m(t)dt+\sigma(t)\epsilon_t$ with instantaneous drift $m$ and instantaneous standard deviation $\sigma$.
Functionality exists for short rate models include:
- 1-factor Vasicek
- Cox Ingersoll Ross
- 2-factor Vasicek (a.k.a Hull-White)
Here we will use historical inflation data to calibrate a 1-factor short rate model.
The tibble includes monthly CPI and annual inflation back to 1948. For more info see help(inflation)
data(inflation) head(inflation)
inflation %>% ggplot(aes(date, inf_yoy)) + geom_line() + scale_y_continuous(labels = scales::percent) + labs(title = "Historical Inflation", y = "Inflation") + theme_classic()
Now that we have our data, we can calibrate a model to project future inflation. Here I will use the full dataset to fit a Vasicek model, which is mean reverting one-factor short rate model of the form $dr_t = a(b - r_t)dt+\sigma\epsilon_t$ For more see here
inf_parm <- CalVasicek1f(inflation$inf_yoy) inf_parm
The result is a list with the following parameters:
- r0 - the starting value (last value of the data provided)
- a - strength of the mean reversion (annual)
- b - mean reverting level
- v - volatility (annual)
Now that we have our parameters, we can simulate future inflation paths. We will do this with a call to SimVasicek1F.
Let's start by simulating 10 1-year paths. Data is simulated in monthly increments with a default time of 1 year.
inf_sim <- SimVasicek1F(n=10, inf_parm) head(inf_sim) inf_sim %>% bind_rows(tibble(trial = 1:10, time = 0, value = inf_parm$r0)) %>% ggplot(aes(time, value, group = trial)) + geom_line(alpha = 1/2) + scale_y_continuous(labels = scales::percent) + labs(title = "Future Projected Inflation", y = "Inflation") + theme_classic()
If we're interested in projections further out we can adjust the time parameter t.
inf_sim <- SimVasicek1F(n=10, inf_parm, t = 60) inf_sim inf_sim %>% bind_rows(tibble(trial = 1:10, time = 0, value = inf_parm$r0)) %>% ggplot(aes(time, value, group = trial)) + geom_line(alpha = 1/2) + scale_y_continuous(labels = scales::percent) + labs(title = "Future Projected Inflation", y = "Inflation")
Now we have 10 simulated 5 year paths.
Now let's simulate 1000 paths and view output.
Convenience function quibble provides some summary statistics.
We can see that as we move further out in time, our mean level is getting closer to the mean. There's also greater dispersion
inf_sim <- SimVasicek1F(n=1e3, inf_parm) inf_sim %>% mutate(time = round(time, 3)) %>% group_by(time) %>% summarise(quibble(value)) %>% pivot_wider(names_from = time) %>% mutate_if(is.numeric, ~ scales::percent(., accuracy = .01)) %>% flextable::flextable() %>% flextable::add_header_lines("Inflation Percentiles by Projection Month")
We can also visualize the range of distributions using a funnel graph gg_funnel.
This shows historical values alongside ssimulation quantiles.
This gives us a sense not only of the ranges in the projected outcomes, but also how they compare to historical results.
gg_funnel(inf_sim, tail(inflation$inf_yoy, 61), inf_parm$r0, variable_name = "Inflation", time_offset = 2022)
Equities
Supported modeling options for equities include:
- [Independent Lognormal]
- [Regime Switching Lognormal]
Independent Lognormal
Assumes monthly log-returns are independent and normally distribution.
The simulation output includes monthly returns, cumulative returns and accumulated value (wealth ratio)
SimEquityILN()
If we wanted to calibrate the parameters using actual data we could utilize the equity dataset. This includes monthly (continuous) returns back to December 1927.
data(equity) head(equity)
We will use the past 20 years to calibrate the lognormal parameters uses maximum likelihood.
iln_param <- CalILN(tail(equity$logreturn, 240))
Here we are assuming a r scales::percent(iln_param[1], accuracy = 0.1) mean and r scales::percent(iln_param[2], accuracy = 0.1) volatility.
(note: the mean is expressed as the annual arithmetic return.)
We can now use those parameters to simulate.
sim_iln <- SimEquityILN(n=1000, iln_param) sim_iln %>% mutate(time = round(time, 3)) %>% group_by(time) %>% summarise(quibble(cum_return)) %>% pivot_wider(names_from = time) %>% mutate_if(is.numeric, ~ scales::percent(., accuracy = .01)) %>% flextable::flextable() %>% flextable::add_header_lines("Inflation Percentiles by Projection Month")
Regime Switching Lognormal
Historical data shows evidence of volatility "bunching", a feature not captured in the independent lognormal model.
equity %>% filter(date > "1960-01-01") %>% mutate(`Volatility` = zoo::rollapply(logreturn, 12, sd, fill = NA) * sqrt(12)) %>% rename(`Total Return` = logreturn) %>% pivot_longer(c(-date)) %>% ggplot(aes(date, value, group = name)) + geom_line(aes(linetype=name)) + scale_y_continuous(labels = scales::percent) + labs(title = "S&P 500 Monthly Total Returns and Annual Volatility", x = "Year", y = NULL, caption = "volatility = 12-month moving standard deviation of log returns") + theme(legend.title = element_blank(), plot.caption = element_text(hjust = 0)) + guides(linetype = guide_legend(reverse = T))
A regime switching model allows for the modeling of this dynamic using distinct parameters in different states or "regimes".
Within a given state ("regime"), log returns are independent, normally distributed, but the markov process of regime switching allows for volatility clustering.
In addition to regime mean and volatility of returns, probabilities of transitioning between regimes is required.
Below we will simulate a single trial with 1000 monthly timesteps.
set.seed(1) x <- SimEquityRSLN(n = 1, pswitch = c(.011, .059), means = c(.008, -.011), vols = c(.039, .113), .detail = T, t = 1e3) x %>% ggplot(aes(x = seq(1:1e3), y = state)) + geom_line() + scale_y_continuous(breaks = c(0, 1)) + labs(x = "Month")
You can see the shifts between the baseline (0 state) and higher volatility states The unconditional probability of being in state 1 is $\frac{p_{1,2}}{p_{1,2} + p_{2,1}}$ or ~15% under these assumptions
