**unstarched» R**, and kindly contributed to R-bloggers)

Since the last release of rugarch on CRAN (ver 1.0-16), there have been many changes and new features in the development version of the package (ver 1.01-5). First, development of the package (and svn) has been moved to google code from r-forge. Second, the package now features exclusive use of xts based time series for input data and also outputs some of the results as xts as well. The sections that follow highlight some of the key changes to the package which I hope will make it easier to work with.

## Extractor Methods

### sigma, fitted and residuals

The main extractor method is no longer the **as.data.frame**, but instead, **sigma** and **fitted** will now extract the conditional sigma (GARCH) and mean (ARFIMA) values from all objects. The old extractor methods are now mostly deprecated, with the exception of certain classes where it still makes sense to use them (i.e. uGARCHroll, uGARCHdistribution, uGARCHboot).

library(rugarch) data(sp500ret) fit = ugarchfit(ugarchspec(), sp500ret, out.sample = 10) c(is(sigma((fit))), is(fitted(fit)), is(residuals(fit))) ## [1] 'xts' 'xts' 'xts' ## plot(xts([email protected]$modeldata$data, [email protected]$modeldata$index), auto.grid = FALSE, minor.ticks = FALSE, main = 'S&P500 Conditional Mean') lines(fitted(fit), col = 2) grid()

plot(xts(abs([email protected]$modeldata$data), [email protected]$modeldata$index), auto.grid = FALSE, minor.ticks = FALSE, main = 'S&P500 Conditional Sigma', col = 'grey') lines(sigma(fit), col = 'steelblue') grid()

Apart from getting the nice xts charts, **rugarch** can now handle any type of time-formatted data (which can be coerced to xts) including intraday.

A key change has also been made to the output of the forecast class (uGARCHforecast) which merits particular attention. Because **rugarch** allows to combine both rolling and unconditional forecasts, this creates a rather challenging problem in how to meaningfully output the results, with simple extractor methods such as **sigma** and **fitted**. The output in this case will be an n.ahead by (n.roll+1) matrix, where the column headings are now the *T+0* dates, since they will always exist to use within a forecast framework, whilst the row names are labelled as *T+1*, *T+2*, …,*T+n.ahead*. The following example illustrates, and shows some simple solutions to deal with portraying n.ahead dates which are completely in the future (i.e. not in the available out of sample testing period).

f = ugarchforecast(fit, n.ahead = 25, n.roll = 9) sf = sigma(f) head(sf) ## 2009-01-15 2009-01-16 2009-01-20 2009-01-21 2009-01-22 2009-01-23 2009-01-26 2009-01-27 2009-01-28 2009-01-29 ## T+1 0.02176 0.02078 0.02587 0.02730 0.02646 0.02519 0.02400 0.02301 0.02388 0.02487 ## T+2 0.02171 0.02073 0.02580 0.02722 0.02638 0.02512 0.02393 0.02295 0.02382 0.02480 ## T+3 0.02166 0.02069 0.02573 0.02714 0.02630 0.02505 0.02387 0.02289 0.02375 0.02473 ## T+4 0.02160 0.02064 0.02565 0.02706 0.02623 0.02498 0.02380 0.02283 0.02369 0.02466 ## T+5 0.02155 0.02059 0.02558 0.02697 0.02615 0.02491 0.02374 0.02277 0.02363 0.02459 ## T+6 0.02150 0.02054 0.02550 0.02689 0.02607 0.02484 0.02368 0.02271 0.02356 0.02452

Therefore, the column headings are the *T+0* dates i.e. the date at which the forecast was made. There are 10 columns since the n.roll option is zero based. To create an xts representation of only the rolling output is very simple:

print(sf[1, ]) ## 2009-01-15 2009-01-16 2009-01-20 2009-01-21 2009-01-22 2009-01-23 2009-01-26 2009-01-27 2009-01-28 2009-01-29 ## 0.02176 0.02078 0.02587 0.02730 0.02646 0.02519 0.02400 0.02301 0.02388 0.02487 ## print(f.roll < - as.xts(sf[1, ])) ## [,1] ## 2009-01-15 0.02176 ## 2009-01-16 0.02078 ## 2009-01-20 0.02587 ## 2009-01-21 0.02730 ## 2009-01-22 0.02646 ## 2009-01-23 0.02519 ## 2009-01-26 0.02400 ## 2009-01-27 0.02301 ## 2009-01-28 0.02388 ## 2009-01-29 0.02487

This gives the *T+0* dates. If you want to generate the actual **T+1** rolling dates, use the **move** function on the *T+0* dates which effectively moves the dates 1-period ahead, and appends an extra 1 day to the end (which is not in the sample date set, therefore this is a ‘generated’ date):

print(f.roll1 < - xts(sf[1, ], move(as.POSIXct(colnames(sf)), by = 1))) ## [,1] ## 2009-01-16 0.02176 ## 2009-01-20 0.02078 ## 2009-01-21 0.02587 ## 2009-01-22 0.02730 ## 2009-01-23 0.02646 ## 2009-01-26 0.02519 ## 2009-01-27 0.02400 ## 2009-01-28 0.02301 ## 2009-01-29 0.02388 ## 2009-01-30 0.02487

For the n.ahead forecasts, I have included in the **rugarch** package another simple date function called **generatefwd** which can be used to generate future non-weekend dates. To use this, you will need to know that the model slot of any class object in **rugarch** which took a method with a data option will keep that data, its format and periodicity, which can be used as follows:

args(generatefwd) ## function (T0, length.out = 1, by = 'days') ## NULL print([email protected]$modeldata$period) ## Time difference of 1 days ## i = 1 DT0 = generatefwd(T0 = as.POSIXct(colnames(sf)[i]), length.out = 25, by = [email protected]$modeldata$period) print(fT0 < - xts(sf[, i], DT0)) ## [,1] ## 2009-01-16 0.02176 ## 2009-01-19 0.02171 ## 2009-01-20 0.02166 ## 2009-01-21 0.02160 ## 2009-01-22 0.02155 ## 2009-01-23 0.02150 ## 2009-01-26 0.02145 ## 2009-01-27 0.02139 ## 2009-01-28 0.02134 ## 2009-01-29 0.02129 ## 2009-01-30 0.02124 ## 2009-02-02 0.02119 ## 2009-02-03 0.02114 ## 2009-02-04 0.02109 ## 2009-02-05 0.02104 ## 2009-02-06 0.02099 ## 2009-02-09 0.02094 ## 2009-02-10 0.02089 ## 2009-02-11 0.02084 ## 2009-02-12 0.02079 ## 2009-02-13 0.02075 ## 2009-02-16 0.02070 ## 2009-02-17 0.02065 ## 2009-02-18 0.02060 ## 2009-02-19 0.02055

The same applies to the fitted method. This concludes the part on the uGARCHforecast class and the changes which have been made. More details can always be found in the documentation (if you are wondering how to obtain documentation for an S4 class object such as uGARCHforecast, you need to append a ‘-class’ to the end and enclose the whole thing in quotes e.g. help(‘uGARCHforecast-class’)).

### quantile and pit

The **quantile** method extracts the conditional quantiles of a **rugarch** object, subject to an extra option (probs) as in the S3 class method in *stats*, while **pit** calculates and returns the probability integral transformation of a fitted, filtered or rolling object (objects guaranteed to have ‘realized’ data to work with).

head(quantile(fit, c(0.01, 0.025, 0.05))) ## q[0.01] q[0.025] q[0.05] ## 1987-03-10 -0.02711 -0.02276 -0.01901 ## 1987-03-11 -0.02673 -0.02247 -0.01881 ## 1987-03-12 -0.02549 -0.02141 -0.01791 ## 1987-03-13 -0.02449 -0.02058 -0.01722 ## 1987-03-16 -0.02350 -0.01972 -0.01647 ## 1987-03-17 -0.02271 -0.01903 -0.01586 ## head(pit(fit)) ## pit ## 1987-03-10 0.7581 ## 1987-03-11 0.4251 ## 1987-03-12 0.5973 ## 1987-03-13 0.3227 ## 1987-03-16 0.2729 ## 1987-03-17 0.9175

These should prove particularly useful in functions which require the quantiles or PIT transformation, such as the risk (**VaRTest** and **VaRDurTest**) and misspecification tests (**BerkowitzTest** and **HLTest**).

## Distribution Functions

rugarch exports a set of functions for working with certain measures on the conditional distributions included in the package. Some of these functions have recently been re-written to take advantage of vectorization, whilst others re-written to take advantage of analytical representations (e.g. as regards to skewness and kurtosis measures for the sstd and jsu distributions). In the latest release, I’ve also included some functions which provide for a graphical visualization of the different distributions with regards to their higher moment features, and demonstrated below:

distplot('nig') ## Warning: skew lower bound below admissible region...adjusting to ## distribution lower bound.

The **distplot** function, available for all skewed and/or shaped distributions provides a visual 3D representation of the interaction of the skew and shape parameters in determining the Skewness and Kurtosis. In cases where the distribution is only skewed or shaped, a simpler 2D plot is created as in the case of the *std* distribution:

distplot('std')

Another interesting plot is that of a distribution’s ‘authorized domain’, which shows the region of Skewness-Kurtosis for which a density exists. This is related to the Hamburger moment problem, and the maximum attainable Skewness (S) given kurtosis (K) ( see Widder (1946) ) . The **skdomain** function returns the values needed to visualize these relationships, and demonstrated below:

# plot the first one XYnig = skdomain('nig', legend = FALSE) XYsstd = skdomain('sstd', plot = FALSE) XYhyp = skdomain('ghyp', plot = FALSE, lambda = 1) XYjsu = skdomain('jsu', plot = FALSE) # the values returned are the bottom half of the domain (which is # symmetric) lines(XYjsu$Kurtosis, XYjsu$Skewness, col = 3, lty = 2) lines(XYjsu$Kurtosis, -XYjsu$Skewness, col = 3, lty = 2) lines(XYsstd$Kurtosis, XYsstd$Skewness, col = 4, lty = 3) lines(XYsstd$Kurtosis, -XYsstd$Skewness, col = 4, lty = 3) lines(XYhyp$Kurtosis, XYhyp$Skewness, col = 5, lty = 4) lines(XYhyp$Kurtosis, -XYhyp$Skewness, col = 5, lty = 4) legend('topleft', c('MAX', 'NIG', 'JSU', 'SSTD', 'HYP'), col = c(2, 'steelblue', 3, 4, 5), lty = c(1, 1, 2, 3, 4), bty = 'n')

From the plot, it is clear that the skew-student has the widest possible combination of skewness and kurtosis for values of kurtosis less than ~6, whilst the NIG has the widest combination for values greater than ~8. It is interesting to note that the Hyperbolic distribution is not defined for values of kurtosis greater than ~9.

## Model Reduction

The **reduce** function eliminates non-significant coefficients (subject to a pvalue cut-off argument) from a model by fixing them to zero (in **rugarch** this is equivalent to eliminating them) and re-estimating the model with the remaining parameters as the following example illustrates:

spec = ugarchspec(mean.model = list(armaOrder = c(4, 4)), variance.model = list(garchOrder = c(2, 2))) fit = ugarchfit(spec, sp500ret[1:1000, ]) round([email protected]$robust.matcoef, 4) ## Estimate Std. Error t value Pr(>|t|) ## mu 0.0005 0.0004 1.250 0.2111 ## ar1 0.2250 0.0164 13.677 0.0000 ## ar2 1.6474 0.0175 94.046 0.0000 ## ar3 -0.0661 0.0188 -3.517 0.0004 ## ar4 -0.8323 0.0164 -50.636 0.0000 ## ma1 -0.2328 0.0065 -35.825 0.0000 ## ma2 -1.6836 0.0039 -430.423 0.0000 ## ma3 0.0980 0.0098 9.953 0.0000 ## ma4 0.8426 0.0075 112.333 0.0000 ## omega 0.0000 0.0000 1.527 0.1267 ## alpha1 0.2153 0.1265 1.702 0.0888 ## alpha2 0.0000 0.0000 0.000 1.0000 ## beta1 0.3107 0.1149 2.705 0.0068 ## beta2 0.3954 0.0819 4.828 0.0000 ## # eliminate alpha2 fit = reduce(fit, pvalue = 0.1) print(fit) ## ## *---------------------------------* ## * GARCH Model Fit * ## *---------------------------------* ## ## Conditional Variance Dynamics ## ----------------------------------- ## GARCH Model : sGARCH(2,2) ## Mean Model : ARFIMA(4,0,4) ## Distribution : norm ## ## Optimal Parameters ## ------------------------------------ ## Estimate Std. Error t value Pr(>|t|) ## mu 0.00000 NA NA NA ## ar1 0.23153 0.006583 35.1698 0.000000 ## ar2 1.63558 0.018294 89.4050 0.000000 ## ar3 -0.07264 0.018277 -3.9743 0.000071 ## ar4 -0.82110 0.014031 -58.5206 0.000000 ## ma1 -0.23883 0.009659 -24.7253 0.000000 ## ma2 -1.68338 0.005401 -311.6821 0.000000 ## ma3 0.11434 0.017035 6.7121 0.000000 ## ma4 0.83198 0.013140 63.3179 0.000000 ## omega 0.00000 NA NA NA ## alpha1 0.14532 0.008997 16.1530 0.000000 ## alpha2 0.00000 NA NA NA ## beta1 0.26406 0.086397 3.0564 0.002240 ## beta2 0.58961 0.088875 6.6342 0.000000 ## ## Robust Standard Errors: ## Estimate Std. Error t value Pr(>|t|) ## mu 0.00000 NA NA NA ## ar1 0.23153 0.061687 3.7533 0.000175 ## ar2 1.63558 0.093553 17.4829 0.000000 ## ar3 -0.07264 0.024012 -3.0251 0.002485 ## ar4 -0.82110 0.028395 -28.9168 0.000000 ## ma1 -0.23883 0.037806 -6.3172 0.000000 ## ma2 -1.68338 0.021478 -78.3771 0.000000 ## ma3 0.11434 0.066989 1.7068 0.087851 ## ma4 0.83198 0.052093 15.9711 0.000000 ## omega 0.00000 NA NA NA ## alpha1 0.14532 0.057150 2.5428 0.010996 ## alpha2 0.00000 NA NA NA ## beta1 0.26406 0.129725 2.0356 0.041793 ## beta2 0.58961 0.103879 5.6759 0.000000 ## ## LogLikelihood : 3091 ## #####

## References

Hansen, P. R., Lunde, A., & Nason, J. M. (2011). The model confidence set. Econometrica, 79(2), 453-497.

Widder, D. V. (1959). The Laplace Transform. 1946. Princeton University Press, Princeton, New Jersey. Supplementary Bibliography Is. McShane, EJ,“ A canonical form for antiderivatives,” Illinois Jour, of Math, 3, 334-351.

**leave a comment**for the author, please follow the link and comment on their blog:

**unstarched» R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...