Some R Code Examples

This is just some miscellaneous code relating to implied volatility computations, and GARCH modeling.
 library(fBasics)
 library(fGarch)
 library(RQuantLib)
 setwd("c:/Work/Lectures/Courses/HIT_Shenzhen")
 source("financetools_funs.R")

############# Example for VIX implied volatility computations.
###### First just look at VIX daily
 VIX_dc <- get.stock.price(
         "^VIX",start.date=c(1,1,2008),stop.date=c(12,31,2014))
plot(VIX_dc,type='l',
   xlab="Jan 1, 2008 to Dec 31, 2014",ylab="daily closes",
   main="VIX Daily Prices")
   
###### To compute implied volatility, use EuropeanOptionImpliedVolatility
help(EuropeanOptionImpliedVolatility)

## This yields:
## Usage
## EuropeanOptionImpliedVolatility(type, value,
##	underlying, strike, dividendYield, riskFreeRate, maturity, volatility)

##type - A string with one of the values call or put
##value - Value of the option (used only for ImpliedVolatility calculation)
##underlying - Current price of the underlying stock
##strike - Strike price of the option
##dividendYield - Continuous dividend yield (as a fraction) of the stock
##riskFreeRate - Risk-free rate
##maturity - Time to maturity (in fractional years)
##volatility - Initial guess for the volatility of the underlying stock
 
###### Now, get prices from finance.yahoo.com   
## On this day, the VIX was 15.29 (That's the price of the underlying.)
## For the June 17 expiry, and a strike of 15, 
## the bid was 0.95 and the ask was 1.05

## The dividend rate is 0
## The risk-free (now) is very low (ZIRP): 0.005
## The time to maturity is
  Tmat <- as.numeric(difftime("2014-06-17","2014-06-09"))/365.24
## We need a guess for starting; use 0.6 or 0.7 (VIX has large volatility.
## Now call EuropeanOptionImpliedVolatility


###### Examples of GARCH modeling
###### IBM monthly returns for GARCH examples
 IBM_mc <- get.stock.price(
         "IBM",start.date=c(1,1,2009),stop.date=c(12,31,2013),freq="m")
 IBM_mr <- log.ratio(IBM_mc)
plot(IBM_mr,type='l',
   xlab="Jan 1, 2009 to Dec 31, 2013",ylab="monthly log returns",
   main="IBM Monthly Returns")
abline(0,0)
plot(acf(IBM_mr,plot=FALSE)[2:30],main="ACF for IBM Monthly Returns")
Box.test(IBM_mr,lag=5,type="Ljung")

## Box test p-value is 0.34; decide not to fit a conditional mean 
## (i.e., ARMA(0,0) is OK)
#################  GARCH no mean fitting
garchFit(data=IBM_mr, formula=~garch(1,1), trace=FALSE)

## This produces p-values that indicate that this model is not good.
## Unfortunately, we conclude that a ARMA(p,q)+GARCH(m,s) model just does
## not fit the IBM returns.

########  Now let's try another series; this time of an index.
####### S&P 500 daily
GSPC_dc <- get.stock.price(
         "^GSPC",start.date=c(1,1,2008),stop.date=c(12,31,2013))
GSPC_dr <- log.ratio(GSPC_dc)
plot(GSPC_dr,type='l',
   xlab="Jan 1, 2008 to Dec 31, 2013",ylab="daily log returns",
   main="S&P 500 Daily Returns")

plot(acf(GSPC_dr,plot=FALSE)[2:30],main="ACF for S&P 500 Daily Returns")

Box.test(GSPC_dr,lag=5,type="Ljung")

pacf(GSPC_dr)

## This yields highly significant autocorrelations and partial autocorrelations
## out beyond lag ten.  Let's just try a conservative model; p=1 and q=1.
## Fit ARMA(1,1) for mean:

armafit<-arima(GSPC_dr,order=c(1,0,1))
ressq<-armafit$residuals^2
plot(acf(ressq,plot=FALSE)[2:30],main="ACF for Squared Residuals")

Box.test(ressq,lag=5,type="Ljung")

pacf(ressq)

## This yields highly significant autocorrelations and partial autocorrelations
## out beyond lag ten.  Let's just try a conservative model; m=1 and s=1.
## Fit GARCH(1,1) for variance of the ARMA(1,1) residuals

garchFit(data=GSPC_dr, formula=~arma(1,1)+garch(1,1), trace=FALSE)

## This yields 
##         mu          ar1          ma1        omega       alpha1        beta1  
## 2.5664e-04   6.7521e-01  -7.3781e-01   2.2077e-06   1.1017e-01   8.7843e-01  
## All highly significant except for mu.
## The normalized log likelihood is 3.054381.

## Maybe we should try other models.  Try GARCH(2,2) for variance of the 
## ARMA(1,1) residuals.

garchFit(data=GSPC_dr, formula=~arma(1,1)+garch(2,2), trace=FALSE)

## This yields 
## Coefficient(s):
##         mu          ar1          ma1        omega       alpha1       alpha2  
## 2.6230e-04   6.6680e-01  -7.3516e-01   3.5089e-06   1.0000e-08   1.5194e-01  
##      beta1        beta2  
## 8.1873e-01   9.9953e-03 
## ar1, ma1, and alpha2 are highly significant.
## The normalized log likelihood is 3.064624.

garchFit(data=GSPC_dr, formula=~arma(2,2)+garch(2,2), trace=FALSE)

## This yields 
## Coefficient(s):
##          mu          ar1          ar2          ma1          ma2        omega  
##  2.6339e-04   1.5928e-01   5.1659e-01  -2.3290e-01  -5.2597e-01   3.4684e-06  
##      alpha1       alpha2        beta1        beta2  
##  1.0000e-08   1.5041e-01   8.3064e-01   1.0000e-08 
## ar2, ma2, and alpha2 are highly significant.
## The normalized log likelihood is 3.066124.

## There is not much difference in the ARMA(1,1)+ GARCH(2,2) model and the
## ARMA(2,2)+ GARCH(2,2) model.
## I may choose the ARMA(2,2)+ GARCH(2,2) model, however.  Its normalized
## log likelihood is slightly larger.

####### Now try same thing for S&P 500 monthly
GSPC_mc <- get.stock.price(
         "^GSPC",start.date=c(1,1,2008),stop.date=c(12,31,2013),freq="m")
GSPC_mr <- log.ratio(GSPC_mc)
plot(GSPC_mr,type='l',
   xlab="Jan 1, 2008 to Dec 31, 2013",ylab="monthly log returns",
   main="S&P 500 Monthly Returns")

plot(acf(GSPC_mr,plot=FALSE)[2:30],main="ACF for S&P 500 Monthly Returns")

Box.test(GSPC_mr,lag=5,type="Ljung")

###  Decide must fit ARMA(1,1) for mean

garchFit(data=GSPC_mr, formula=~arma(1,1)+garch(1,1), trace=FALSE)