Here you will find daily news and tutorials about R, contributed by over 750 bloggers.
There are many ways to follow us - By e-mail:On Facebook: If you are an R blogger yourself you are invited to add your own R content feed to this site (Non-English R bloggers should add themselves- here)

Adam Duncan, December 2012 Also avilable on R-bloggers.com

Prelude

These posts are written with dual purpose: 1) Hopefully provide some insight or inspiration into a topical issue in finance from a practioners perspective, and 2) show how to use R to craft an analysis and produce nice output. The posts are written in a “walkthrough” style. All of the source files are available at my github repository. Feel free to use the files for your own projects. I only ask that you credit appropriately.

Introduction: How big is a financial “crisis”?tap_a=5644-dce66f&tap_s=10907-287229"SP500","DJIA","DEXUSEU","DSWP5","DCOILWTICO","DEXMXUS")for(i in1:length(tickers)){
suppressWarnings(getSymbols(tickers[i], src ="FRED", auto.assign =TRUE))}

Let’s just focus on one data set for the moment. For example, let’s choose USDMXN. For this post, we’ll just construct the analysis for USDMXN. In the next post, we’ll run through many variables and try to draw some more conclusions about the magnitude of past crises.

data <- DEXMXUS # Be careful when selecting data. For example, USDMXN goes UP in a crisis...
data <-1/data # Since we're using USDMXN as an example will flip it to MXNUSD so it goes down in a crisis.

Next we’ll need to calculate a return series for the data object…

ret <- periodReturn(data, period ="daily", subset =NULL, type ="log")

Similarly, we’ll need to calculate a rolling mean, std.dev., skewness and kurtosis for our series, given our VaR lookback length choice: var.yrs.:

ret.mu <- rollapply(ret, width = hist.len, FUN = mean)
ret.sig <- rollapply(ret, width = hist.len, FUN =function(ret) apply(X = ret,
MARGIN =2, FUN = sd))
ret.skew <- rollapply(ret, width = hist.len, FUN =function(ret) skewness(ret,
method ="moment"))
ret.kurt <- rollapply(ret, width = hist.len, FUN =function(ret) kurtosis(ret,
method ="moment"))

We can now use rollapply to calculate the daily VaR we would have seen each day in the sample data. I choose a historical VaR method to keep things simple and a single underlying. You can easily choose “gaussian” and pass in skewness and kurtosis, if you want to make it more consistent with how practitioners make parametric VaR estimates. Since we are using VaR as a measuring stick to measure the intesnity of crises, parametric VaR isn’t completely necessary. Having said that, it will affect the ratio of actual realized returns to a priori VaR estimates. With adjustments for skewness and kurtosis, ratios will be smaller, all else equal. That is, crises events will seem less severe when the distribution is adjusted for skewness and excess kurtosis.

Ok, let’s examine the results of all that to make sure things are looking like they should. I always like to pause to examine some graphics along the way to make sure things are behaving as expected. It’s better to find our early that we have a problem. Some plots:

par(mfrow = c(2,2))
plot.xts(ret.mu, las =1, main ="Rolling Mean Returns", cex.main =0.9)
plot.xts(ret.sig * sqrt(biz.days)*100, las =1, main ="Rolling Vol", cex.main =0.9)
plot.xts(ret.skew, las =1, main ="Rolling Skewness", cex.main =0.9)
plot.xts(ret.kurt, las =1, main ="Rolling Kurtosis", cex.main =0.9)

Great. All seems to be in order. One more to go…the rolling VaR estimates.

We better have a look at that handiwork before we move on…

par(mfrow = c(1,1))
plot.xts(ret.var, las =1, main ="Rolling VaR Estimate", cex.main =0.9)

Now that we have the rolling VaR series, we can calculate the ratios of actual returns to VaR estimates. This will give us the magnitude of outcomes in multiples of daily value at risk. We need to trim off the first part of the series that was used to estimate the mean/sd etc. this makes the return subset the same length as ret.var.

first.date <- index(first(ret.var))
ret.sub <- window(ret, from = first.date)

Now, divide the actual returns by the prior day’s VaR etimate. This ratio represents the multiple that the actual negative return was of the VaR estimate. A ratio of -4x means that the realized negative return was 4x larger in magnitude than than the VaR estimate.

Index daily.returns
Min. :1998-11-24 Min. :-11.094
1st Qu.:2002-06-14 1st Qu.: -0.625
Median :2005-12-23 Median : -0.339
Mean :2006-01-16 Mean : -0.486
3rd Qu.:2009-07-29 3rd Qu.: -0.156
Max. :2013-05-01 Max. : -0.001

plot(var.breaks, las =1, main ="Actual Negative Return(t) / Daily VaR Estimate(t-1)",
cex.main =0.9)# line plot of the breaks.

For USDMXN, we can see that 2008 was 11x daily value at risk. Other, smaller crisis episodes occur at around 4-5x daily VaR. Because we are using 5 years of history for the VaR calculation, we miss the de-pegging of USDMXN that happened in 1994. That was on the order of 25x VaR.

Here is a histogram of the VaR breaks for USDMXN:

hist(var.breaks, main ="Histogram of VaR Breaks (Multiples of Daily VaR", cex.main =0.9)# histogram of the VaR breaks.

Now we can just loop through each of the data series we imported and run the same analysis for each one. Then we’ll examine the individual results as well as aggregate them so we can draw some conclusions about how big, in VaR terms, past crises have been.

A Couple Workhorse Functions

First, here is the workhorse function for the analysis. This is the function I use most when I want to examine VaR breaks for a given time series. This function takes raw time series data in xts format (x), and returns one of a number of items depending on the value of ‘rtype’ passed to the function:

‘rtype’ can take on the following valueS: “vars”,”returns”,”ratios”,”vol”,”skew”,”kurt”,”mean” ‘window’ is the length of the lookback period in days. ’p’ is the quantile desired (e.g. .95, .99, etc.) ‘inv.var’ sets whether you want positive or negative VaRs returned from VaR. See Performance Analytics::VaR Based on user input for rtype, the function takes the appropriate action. I’m only going to support the “historical” VaR method for now.

varFun<-function(x,p,window,method,inv.var,rtype){
ret<-periodReturn(x,period="daily",subset=NULL,type="log")if(rtype=="vars"){
ret.var<-rollapply(ret,width=window,
FUN=function(ret) VaR(R=ret,
p=p,
method="historical",
clean="none",
portfolio_method="single",
weights=NULL,#mu=ret.mu,#sigma=ret.sig,#m3=ret.skew,#skewness of series#m4=ret.kurt,#kurtosis of series
invert=inv.var),
align="right")
ret.var<-as.xts(ret.var)
names(ret.var)<-c("VaR")return(ret.var)}else{if(rtype=="returns"){
ret<-as.xts(ret)
names(ret)<-c("returns")return(ret)}else{if(rtype=="ratios"){
ret.var<-rollapply(ret,width=window,
FUN=function(ret) VaR(R=ret,
p=p,
method="historical",
clean="none",
portfolio_method="single",
weights=NULL,#mu=ret.mu,#sigma=ret.sig,#m3=ret.skew,#skewness of series#m4=ret.kurt,#kurtosis of series
invert=inv.var),
align="right")
first.date<-index(first(ret.var))
ret.sub<-window(ret,from=first.date)
act.var.ratio<-as.xts(ret.sub/ret.var)
names(act.var.ratio)<-c("ratios")return(act.var.ratio)}else{if(rtype=="vol"){
ret.sig<-rollapply(ret,width=window,FUN=function(ret) apply(X=ret,MARGIN=2, FUN=sd))
names(ret.sig)<-c("vols")return(ret.sig)}else{if(rtype=="skew"){
ret.skew<-as.xts(rollapply(ret,width=window,FUN=function(ret) skewness(ret,method="moment")))
names(ret.skew)<-c("skew")return(ret.skew)}else{if(rtype=="kurt"){
ret.kurt<-as.xts(rollapply(ret,width=window,FUN=function(ret) kurtosis(ret,method="moment")))
names(ret.kurt)<-c("kurtosis")return(ret.kurt)}else{if(rtype=="mean"){
ret.mu<-as.xts(rollapply(ret,width=window,FUN=mean))
names(ret.mu)<-c("means")return(as.xts(ret.mu))}else{return(Null)}}}}}}}}

We can now assign some variable, v, to hold the VaR breaks for the variable “data”:

v <- varFun(data,0.95,1260,"historical",FALSE,"ratios")

And then plot them…

plot.xts(subset(v, subset =(v <0)), las =1, main ="Rolling Historical VaR Breaks",
ylab ="Raito of Actual (Neg) Returns to VaR")

If you just have a bunch of tickers and want to see a summary of the VaR breaks, this will do it. This function loops through each of the data elements in ‘tickers’ and calculates some summary information. It returns a data frame of the results.

Note: This is the VERY slow way to do this. It’s better to import the data first, and then past a list of raw data variables to this function. Getting the data on the fly is slow and then you don’t have use of the variables on exit. I’ll update this summary function in the next post.

vsum<-function(tickers,p,window,method,inv.var){
n<-length(tickers)
out<-data.frame()for(i in1:n){
data<-getSymbols(tickers[i],src="FRED",auto.assign=FALSE)
ret<-periodReturn(data,period="daily",subset=NULL,type="log")
ret.mu<-rollapply(ret,width=hist.len,FUN=mean)
ret.sig<-rollapply(ret,width=hist.len,FUN=function(ret) apply(X=ret,MARGIN=2, FUN=sd))
ret.skew<-rollapply(ret,width=hist.len,FUN=function(ret) skewness(ret,method="moment"))
ret.kurt<-rollapply(ret,width=hist.len,FUN=function(ret) kurtosis(ret,method="moment"))
ret.var<-rollapply(ret,width=window,
FUN=function(ret) VaR(R=ret,
p=p,
method=method,
clean="none",
portfolio_method="single",
weights=NULL,
mu=ret.mu,
sigma=ret.sig,
m3=ret.skew,#skewness of series
m4=ret.kurt,#kurtosis of series
invert=inv.var),
align="right")
first.date<-index(first(ret.var))
ret.sub<-window(ret,from=first.date)
act.var.ratio<-ret.sub/ret.var
out[i,1]<-tickers[i]
out[i,2]<-max(act.var.ratio)
out[i,3]<-mean(act.var.ratio)
out[i,4]<-quantile(act.var.ratio,p,na.rm=TRUE)}
names(out)<-c("name","max.ratio","mean.ratio","Q.p")return(out)}

Next Up: Running this analysis across many assets

Related

To leave a comment for the author, please follow the link and comment on their blog: From Guinness to GARCH.