Introduction to R for Quantitative Finance: R is a statistical computing language that's ideal for answering quantitative finance questions. This book gives you both theory and practice, all in clear language with stacks of real-world examples. Ideal for R beginners or expert alike.

We use the auto.arima function provided by the forecast package to identify the optimal model and estimate the coefficients in one step. The function takes several arguments besides the return series (hp_ret). By specifying stationary = TRUE,we restrict the search to stationary models. In a similar vein, seasonal = FALSE restricts the search to non-seasonal models. Furthermore, we select the Akaike information criteria as the measure of relative quality to be used in model selection.

> mod <- auto.arima(hp_ret, stationary = TRUE, seasonal = FALSE,+   ic="aic")

To determine the fitted coefficient values, we query the model output.

> mod
Series: hp_ret
ARIMA(2,0,0) with non-zero mean 

Coefficients:
         ar1     ar2  intercept
      0.2299  0.3491     0.4345
s.e.  0.0573  0.0575     0.1519

sigma^2 estimated as 1.105:  log likelihood=-390.97
AIC=789.94   AICc=790.1   BIC=804.28

An AR(2) process seems to fit the data best, according to Akaike's Information Criteria. For visual confirmation, we can plot the partial autocorrelation function using the command pacf. It shows non-zero partial autocorrelations until lag two, hence an AR process of order two seems to be appropriate. The two AR coefficients, the intercept (which is actually the mean if the model contains an AR term), and the respective standard errors are given. In the following example, they are all significant at the 5% level since the respective confidence intervals do not contain zero:

> confint(mod)
              2.5 %    97.5 %
ar1       0.1174881 0.3422486
ar2       0.2364347 0.4617421
intercept 0.1368785 0.7321623

If the model contains coefficients that are insignificant, we can estimate the model anew using the arima function with the fixed argument, which takes as input a vector of elements 0 and NA. NA indicates that the respective coefficient shall be estimated and 0 indicates that the respective coefficient should be set to zero.

A quick way to validate the model is to plot time-series diagnostics using the following command:

> tsdiag(mod)

The output of the preceding command is shown in the following figure:

Our model looks good since the standardized residuals don't show volatility clusters, no significant autocorrelations between the residuals according to the ACF plot, and the Ljung-Box test for autocorrelation shows high p-values, so the null hypothesis of independent residuals cannot be rejected.

To assess how well the model represents the data in the sample, we can plot the raw monthly returns (the thin black solid line) versus the fitted values (the thick red dotted line).

> plot(mod$x, lty = 1, main = "UK house prices: raw data vs. fitted+   values", ylab = "Return in percent", xlab = "Date")
> lines(fitted(mod), lty = 2,lwd = 2, col = "red")

The output is shown in the following figure:

Furthermore, we can calculate common measures of accuracy.

> accuracy(mod)
ME      RMSE      MAE          MPE   MAPE    MASE
0.00120 1.0514    0.8059       -Inf  Inf     0.792980241

This command returns the mean error, root mean squared error, mean absolute error, mean percentage error, mean absolute percentage error, and mean absolute scaled error.

To predict the monthly returns for the next three months (April to June 2013), use the following command:

> predict(mod, n.ahead=3)
$pred
           Apr       May       Jun
2013 0.5490544 0.7367277 0.5439708

$se
          Apr      May      Jun
2013 1.051422 1.078842 1.158658

So we expect a slight increase in the average home prices over the next three months, but with a high standard error of around 1%. To plot the forecast with standard errors, we can use the following command:

> plot(forecast(mod))

A plot of the returns indicates that ARCH effects might exist in the monthly return data.

> plot(intc, main = "Monthly returns of Intel Corporation",+   xlab = "Date", ylab = "Return in percent")

The output of the preceding commands is as shown in the following figure:

We can use statistical hypothesis tests to verify our inkling. Two commonly used tests are as follows:

First, we perform the Ljung-Box test on the first 12 lags of the squared returns using the following command:

> Box.test(coredata(intc^2), type = "Ljung-Box", lag = 12)

   Box-Ljung test

data:  coredata(intc^2)
X-squared = 79.3451, df = 12, p-value = 5.502e-12

We can reject the null hypothesis of no autocorrelations in the squared returns at the 1% significance level. Alternatively, we could employ the LM test from the FinTS package, which gives the same result.

> install.packages("FinTS")
> library("FinTS")
> ArchTest(coredata(intc))

   ARCH LM-test; Null hypothesis: no ARCH effects

data:  coredata(intc)
Chi-squared = 59.3647, df = 12, p-value = 2.946e-08

Both tests confirm that ARCH effects exist in the monthly Intel returns; hence, an ARCH or GARCH model should be employed in modeling the return time series.

The most commonly used GARCH model, and one that is usually appropriate for financial time series as well, is a GARCH(1,1) model. We use the functions provided by the rugarch library for model specification, parameter estimation, backtesting, and forecasting. If you haven't installed the package, use the following command:

> install.packages("rugarch")

Afterwards, we can load the library using the following command:

> library("rugarch")

First, we need to specify the model using the function ugarchspec. For a GARCH(1,1) model, we need to set the garchOrder to c(1,1) and the model for the mean (mean.model) should be a white noise process and hence equal to armaOrder = c(0,0).

> intc_garch11_spec <- ugarchspec(variance.model = list(+   garchOrder = c(1, 1)),+  mean.model = list(armaOrder = c(0, 0)))

The actual fitting of the coefficients by the method of maximum likelihood is done by the function ugarchfit using the model specification and the return data as inputs.

> intc_garch11_fit <- ugarchfit(spec = intc_garch11_spec,+  data = intc)

For additional arguments, see the Help on ugarchfit. The output of the fitted model (use the command intc_garch11_fit) reveals useful information, such as the values of the optimal parameters, the value of the log-likelihood function, and the information criteria.

A useful test for checking the model performance is to do a historical backtest. In a risk model backtest, we compare the estimated VaR with the actual return over the period. If the return is more negative than the VaR, we have a VaR exceedance. In our case, a VaR exceedance should only occur in 1% of the cases (since we specified a 99% confidence level).

The function ugarchroll performs a historical backtest on the specified GARCH model (here the model is intc_garch11_spec). We specify the backtest as follows:

The following command shows all the preceding points for a backtest:

> intc_garch11_roll <- ugarchroll(intc_garch11_spec, intc,+   n.start = 120, refit.every = 1, refit.window = "moving",+   solver = "hybrid", calculate.VaR = TRUE, VaR.alpha = 0.01,+   keep.coef = TRUE)

We can examine the backtesting report using the report function. By specifying the type argument as VaR, the function executes the unconditional and conditional coverage tests for exceedances. VaR.alpha is the tail probability and conf.level is the confidence level on which the conditional coverage hypothesis test will be based.

> report(intc_garch11_roll, type = "VaR", VaR.alpha = 0.01,+   conf.level = 0.99)
VaR Backtest Report
===========================================
Model:            sGARCH-norm
Backtest Length:   312
Data:            

==========================================
alpha:            1%
Expected Exceed:   3.1
Actual VaR Exceed:	5
Actual %:         1.6%

Unconditional Coverage (Kupiec)
Null-Hypothesis:   Correct Exceedances
LR.uc Statistic:   0.968
LR.uc Critical:      6.635
LR.uc p-value:      0.325
Reject Null:      NO

Conditional Coverage (Christoffersen)
Null-Hypothesis:   Correct Exceedances and
               Independence of Failures
LR.cc Statistic:   1.131
LR.cc Critical:      9.21
LR.cc p-value:      0.568
Reject Null:      O

Kupiec's unconditional coverage compares the number of expected versus actual exceedances given the tail probability of VaR, while the Christoffersen test is a joint test of the unconditional coverage and the independence of the exceedances. In our case, despite the actual five exceedances versus an expectation of three, we can't reject the null hypothesis that the exceedances are correct and independent.

A plot of the backtesting performance can also be generated easily. First, create a zoo object using the extracted forecasted VaR from the ugarchroll object.

> intc_VaR <- zoo(intc_garch11_roll@forecast$VaR[, 1])

We overwrite the index property of the zoo object with rownames (year and month) from this object as well.

> index(intc_VaR) <- as.yearmon(rownames(intc_garch11_roll@forecast$VaR))

We do the same for the actual returns that are also stored in the ugarchroll object.

> intc_actual <- zoo(intc_garch11_roll@forecast$VaR[, 2])
> index(intc_actual) <- as.yearmon(rownames(intc_garch11_roll@forecast$VaR))

Now, we are able to plot the VaR versus the actual returns of Intel using the following commands:

> plot(intc_actual, type = "b", main = "99% 1 Month VaR Backtesting",+   xlab = "Date", ylab = "Return/VaR in percent")
> lines(intc_VaR, col = "red")
> legend("topright", inset=.05, c("Intel return","VaR"), col = c("black","red"), lty = c(1,1))

The following figure shows the output of the preceding command lines:

Now that we can be reasonably sure that our risk model works properly, we can produce VaR forecasts as well. The function ugarchforecast takes as arguments the fitted GARCH model (intc_garch11_fit) and the number of periods for which a forecast should be produced (n.ahead = 12; that is, twelve months).

> intc_garch11_fcst <- ugarchforecast(intc_garch11_fit, n.ahead = 12)

The resulting forecast can be expected by querying the forecast object as shown in the following command lines:

> intc_garch11_fcst
*------------------------------------*
*       GARCH Model Forecast         *
*------------------------------------*
Model: sGARCH
Horizon: 12
Roll Steps: 0
Out of Sample: 0

0-roll forecast [T0=Dec 2008]:
      Series  Sigma
T+1  0.01911 0.1168
T+2  0.01911 0.1172
T+3  0.01911 0.1177
T+4  0.01911 0.1181
T+5  0.01911 0.1184
T+6  0.01911 0.1188
T+7  0.01911 0.1191
T+8  0.01911 0.1194
T+9  0.01911 0.1197
T+10 0.01911 0.1200
T+11 0.01911 0.1202
T+12 0.01911 0.1204

The one-period ahead forecast for the volatility (sigma) is 0.1168. Since we assume a normal distribution, the 99% VaR can be calculated using the 99% quantile (type in qnorm(0.99)) of the standard normal distribution. The one-month 99% VaR estimate for the next period is hence qnorm(0.99)*0.1168 = 0.2717. Hence, with 99% probability the monthly return is above -27%.