107: Stat Introduction to Business and Financial Statistics
Class 21: VaR and Resampling
1
Value at Risk
Simply put, assume you have a $1,000,000
portfolio (very cool! )
You ask: What is the probability that I'd lose
$X (OR MORE) over the next M months?
More precisely, Value at Risk is an estimate of the worst
possible loss an investment could realize over a given time horizon,
under normal market conditions (defined by a given level of
confidence).
To make life simpler, lets ask what will happen over the next 1
month.
For history of VaR see
http://www.riskglossary.com/link/riskmetrics.htm
2
Choice of confidence level 95%
5%
95%
VaR
Investment returns
Normal market conditions the returns that account for 95% of the
distribution of possible outcomes.
Abnormal market conditions the returns that account for the other 5%
of the possible outcomes.
Consider AAPL
5%
4
Example: VaR for AAPL
> getSymbols("AAPL")
[1] "AAPL"
> aapl=monthlyReturn(AAPL)
> length(aapl)
[1] 52
> round(52*.05)
[1] 3
> sort(as.numeric(aapl))
[1] -0.329558190 -0.316639742
[6] -0.078989964 -0.076388889
[11] -0.040694920 -0.027603225
[16] -0.013064272 -0.009097970
[21] 0.022740826 0.029499969
[26] 0.048744570 0.051001821
[31] 0.060530504 0.060722467
[36] 0.079313359 0.079646018
[41] 0.101896439 0.108246678
[46] 0.148470335 0.167215138
[51] 0.214328657 0.237701179
-0.138674598 -0.112900662 -0.088596783
-0.055004859 -0.053404892 -0.050704730
-0.020826845 -0.016124708 -0.013306532
-0.006489744 0.007013780 0.016994875
0.033789621 0.036670416 0.040934811
0.051959325 0.054124356 0.056004687
0.065396230 0.066561812 0.074157787
0.085081920 0.087037647 0.098097152
0.111021277 0.147160008 0.147816349
0.177023850 0.197012938 0.212195122
5
Using the quantile command
Quantiles split the data set into different
percentages.
To find VaR, we want to find the lower 5% of
the returns.
Due to have quantiles are calculated (dont
get me started!) the results wont be the same
as the previous slide.
> quantile(as.numeric(monthlyReturn(AAPL)),.05)
5%
-0.1244989
6
Using the Normal Distribution
Assume that the returns are normally
distirbuted.
Calculate the sample mean and std dev.
Find the lower 5% point on the normal curve.
> qnorm(.05,mean(monthlyReturn(AAPL)),sqrt(var(monthlyReturn(AAPL))))
[1] -0.1509124
7
Using Density Estimation
If we examine a histogram of AAPL returns,
they dont exactly look normal:
8
Density Estimation
Can calculate VaR using the logspline
density estimation command in R.
Compare logspline density estimate versus
normal
normal distribution
> library(logspline)
> fit=logspline(as.numeric(monthlyReturn(AAPL)))
9
In R
> hist(monthlyReturn(AAPL),prob=TRUE)
> lines(x,dnorm(x,m,s))
> lines(x,dlogspline(x,fit))
10
Calculating VaR Using Logspline
We simple need to calculate the logspline
quantile:
> qlogspline(.05,fit)
[1] -0.1668506
11
VaR in Performance Analytics
VaR is a built in function in the
PerformancAnalytics Package
> VaR(monthlyReturn(AAPL),method="historical")
monthly.returns
VaR
-0.1244989
> VaR(monthlyReturn(AAPL),method="gaussian")
monthly.returns
VaR
-0.1491332
Everyone implements this stuff slightly differently so these results
are slightly different than what we got eariler.
12
Confidence Interval for VaR
There is a true VaR but it is unknown.
The VaR we report is just an estimate.
How accurate is it?
We can create a confidence interval for VaR
by resampling.
13
Resampling
The practice of resampling is such an
important and useful concept, that I wanted to
spend another lecture on it.
It does not rely on any modeling
assumptions, and instead, lets the collected
data drive the analysis.
14
Resampling versus Monte Carlo
A Monte Carlo simulation assumes some
inputs to a problem can be simulated as
random variables.
These inputs are then simulated many times,
resulting in many solutions to the problem,
which are then analyzed.
15
MC Example: Flipping a coin
Suppose instead of buying and holding a
stock, you decide to flip a coin each month.
How would this compare to buy and hold?
Clearly for a strongly up-trending stock this is
a stupid idea (ok-its a stupid idea no matter
what).
16
The Code
mikeflip=function(sym) {
stock = getSymbols(sym,auto.assign=FALSE)
sret = monthlyReturn(stock)
buyandhold = prod(1+sret)
svals=1:100
n=length(sret)
for(i in 1:100) {
flips=1.0*(runif(n)>.5)
svals[i]=prod(1+(sret*flips))
}
cat("Coin flip beats buy and hold =
",sum(svals>buyandhold)/100,"\n")
}
17
The Output
> mikeflip("MRK")
Coin flip beats buy
> mikeflip("MRK")
Coin flip beats buy
> mikeflip("AAPL")
Coin flip beats buy
> mikeflip("AAPL")
Coin flip beats buy
> mikeflip("JNJ")
Coin flip beats buy
> mikeflip("JNJ")
Coin flip beats buy
> mikeflip("PRPFX")
Coin flip beats buy
> mikeflip("FMAGX")
Coin flip beats buy
and hold =
0.69
and hold =
0.68
and hold =
0.03
and hold =
0.06
and hold =
0.66
and hold =
0.51
and hold =
0.04
and hold =
0.63
18
Resampling versus Monte Carlo
In resampling, the object of interest is a
statistic which is a function of the data
(x1,x2,x3,..).
This might be
(x1,x2,x3,..) = mean(xs)
(x1,x2,x3,..) = median(xs)
(x1,x2,x3,..) = VaR(xs)
19
How resampling works
One literally resamples the data, that is,
draws random samples of size n from the
observed data set.
For each resampled data set, one can
compute i(x1,x2,x3,..)
This is done say 10000 times, producing
10000 is.
A confidence interval is then calculated.
20
How the CI is calculated
There are several ways to calculate a
confidence interval once 10000 estimates are
obtained:
Empirical CI; sort the values and take the
middle 95%
Normality; mean +/- 1.96*(std dev)
The middle 95% from a density estimate.
21
Example with the mean
The usual 95% confidence interval for a
population mean from a sample of size n is
given by
x t (.025,n-1)
s
n
22
Example: the sample mean
Suppose we want to create a confidence
interval for mean exam score from the
following data:
x t (.025,n-1)
s
n
23
If you didnt know the formula
Resample the returns a large number of
times,
Computing the mean return for each
resample.
resample.
n=10000
vals=1:n
exam=c(78,88,62,90,99,86,85,92,94,75,53,92)
sampsize=length(exam)
for(i in 1:n) {
resampvals=sample(exam,sampsize,replace=TRUE)
vals[i]=mean(resampvals)
}
24
Empirical 95% Confidence Interval
We take the 10000 computed means and sort
them from smallest to largest
The empirical 95% confidence will be the
values between the lower 2.5 percentile and
upper 97.5%
> ord=order(vals)
> svals=vals[ord]
> svals[250]
[1] 74.75
> svals[9750]
[1] 89.66667
Close!
25
CI based on the normal distribution
The mean of the 10000 resampled values +/1.96 * standard deviation
> mean(vals)-1.96*sqrt(var(vals))
[1] 75.44916
> mean(vals)+1.96*sqrt(var(vals))
[1] 90.2535
26
CI based on a density estimate
Fit a density estimate to the resampled
values, then find the inner 95% of that
resulting density.
> fit=logspline(vals)
> qlogspline(.025,fit)
[1] 74.96168
> qlogspline(.975,fit)
[1] 89.67376
27
Your statistical toolbox
Two important techniques in your toolbox
Resampling
Use it whenever you devise a new function of your
data (like Gain-Loss Spread) and what to
understand its variability
Density estimation
Use it to model historical returns, output of a
process, etc.
28
BTW, History of Density Estimation
In statistics, kernel density estimation is a non-parametric way of
estimating the probability density function of a random variable. Kernel
density estimation is a fundamental data smoothing problem where
inferences about the population are made, based on a finite data
sample. In some fields such as signal processing and econometrics it
is also known as the Parzen-Rosenblatt window method, after
Emanuel Parzen and Murray Rosenblatt, who are usually credited with
independently creating it in its current form,[1][2].
29
Recall Value at Risk (Var)
Easiest to look at a picture
30
Finding VaR by Density Estimation
Can calculate VaR using the logspline
density estimation command in R.
> library(logspline)
>
fit=logspline(as.numeric(monthlyReturn(AAPL)))
31
In R
> hist(monthlyReturn(AAPL),prob=TRUE)
> lines(x,dnorm(x,m,s))
> lines(x,dlogspline(x,fit))
32
Calculating VaR Using Logspline
We simply need to calculate the logspline
quantile:
> library(logspline)
>
fit=logspline(as.numeric(monthlyReturn(AAPL)))
> qlogspline(.05,fit)
[1] -0.1668506
33
Pictorially
Area = 5%
-0.1668506
34
Resampling VaR
Code based on density estimation
n=1000
vals=1:n
rets=monthlyReturn(AAPL)
sampsize=length(rets)
fit=logspline(as.numeric(rets))
for(i in 1:n) {
resampvals=rlogspline(sampsize,fit)
vals[i]=VaR(resampvals) ## many ways to compute VaR
}
35
Analyzing the Output
The VaR, and 95% confidence interval
> qlogspline(.05,fit)
[1] -0.1674382
> fit2=logspline(vals)
> -0.2695005
> qlogspline(.025,fit2)
[1] qlogspline(.975,fit2)
[1] -0.07528924
36
Resampling = the Bootstrap
The Bootstrap is another name for
resampling pick yourself up by your
bootstraps
Let the data tell you about itself.
37
Brad Efron
Brad Efron from Stanford invented the bootstrap
(and proved this simple yet elegant idea works).
Efron Dad
Mom Me
38
The Bootstrap in R
Although resampling/bootstrap code is easy
to write, there is a bootstrap routine built into
R.
library(bootstrap)
The function call:
bootstrap(x,nboot,theta,..., func=NULL)
The data
The number
of bootstraps
The function to bootstrap
39
The Mean Example Again
Running this in R
> bfit=bootstrap(exam,10000,mean)
> fit=logspline(bfit$thetastar)
> qlogspline(.025,fit)
[1] 74.90456
> qlogspline(.9755,fit)
[1] 89.5573
See bootstrap package for more ways to compute confidence intervals
40
Bootstrap VaR
Now easy to bootstrap VaR
> getSymbols("AAPL")
[1] "AAPL"
> aaplret=monthlyReturn(AAPL)
> bfit=bootstrap(as.numeric(aaplret),1000,VaR)
> bfits=bfit$thetastar
> fit=logspline(bfits)
> qlogspline(.025,fit)
[1] -0.2552234
> qlogspline(.975,fit)
Note-not a symmetric
[1] -0.07512234
confidence interval-doesnt
> VaR(aaplret)
have to be!
monthly.returns
VaR
-0.1736913
41
Remember the Bootstrap
The bootstrap (resampling) is an incredibly
powerful tool to keep around.
It is easy to explain to people
Enables one to create confidence intervals
for anything you can estimate from the data,
(VaR, Omega, whatever) without worrying
about complex math.
42
Playing with Density Estimation
Risk Measures
We have seen that the standard deviation of
returns is a useful risk measure.
What if we calculate the probability of losing
money
money and use that as a risk measure?
43
Example
Compare AAPL, GOOG and JNJ
By standard deviation
> sqrt(var(aaplret))
monthly.returns
monthly.returns
0.1120790
> sqrt(var(googret))
monthly.returns
monthly.returns
0.1016935
> sqrt(var(jnjret))
monthly.returns
monthly.returns
0.04568969
44
Compare by VaR
Comparing by VaR
> VaR(aaplret)
monthly.returns
VaR
-0.1736913
> VaR(googret)
monthly.returns
VaR
-0.1460762
> VaR(jnjret)
monthly.returns
VaR
-0.0812052
45
Compare by P(ret < 0)
Compute probability of a negative return
> fit.appl=logspline(as.numeric(aaplret))
> fit.goog=logspline(as.numeric(googret))
> fit.jnj=logspline(as.numeric(jnjret))
> plogspline(0,fit.appl)
[1] 0.3064950
> plogspline(0,fit.goog)
[1] 0.4355341
> plogspline(0,fit.jnj)
[1] 0.4694262
Riskiest??
46
Remember Fat Tails
JNJ
AAPL
GOOG
This is why VaR is a useful measure
47
Playing with Density Estimation Random Forecasting
Suppose we took a bunch of historical daily
returns data and fit a density estimate to it.
We then randomly generate 20 (say) returns
from this distribution to predict a stock price
20 days into the future.
How well would it do?
48
The Code
The basic R code
getSymbols("GE",from="2010-01-01",to="2011-03-11")
ret=as.numeric(dailyReturn(GE))
fit=logspline(ret)
samp=rlogspline(20,fit)
newp=last(Cl(GE))*prod(1+samp)
The variable newp is our predicted price 20
days in the future.
49
How does it do?
Here is some R code as a function now
mpredict=function(sym) {
sdata=getSymbols(sym,from="2010-01-01",to="2011-0311",auto.assign=FALSE)
ret=as.numeric(dailyReturn(sdata))
fit=logspline(ret)
samp=rlogspline(20,fit)
newp=last(Cl(sdata))*prod(1+samp)
last20 = Cl(getSymbols(sym,from="2011-03-12",auto.assign=FALSE))
cat("PREDICTED PRICE = ",newp,"ACTUAL PRICE =",last(last20),"\n")
}
50
How does it do?
Here is some output-Ill take that first
prediction
> mpredict("GE")
PREDICTED PRICE =
> mpredict("GE")
PREDICTED PRICE =
> mpredict("GE")
PREDICTED PRICE =
> mpredict("GE")
PREDICTED PRICE =
> mpredict("GE")
PREDICTED PRICE =
> mpredict("GE")
PREDICTED PRICE =
> mpredict("GE")
PREDICTED PRICE =
20.15663 ACTUAL PRICE = 20.19
22.14141 ACTUAL PRICE = 20.19
20.24312 ACTUAL PRICE = 20.19
22.50563 ACTUAL PRICE = 20.19
22.30143 ACTUAL PRICE = 20.19
21.89173 ACTUAL PRICE = 20.19
22.00789 ACTUAL PRICE = 20.19
51
Looks like this idea sucks
Trying SPY
> mpredict("SPY")
PREDICTED PRICE =
> mpredict("SPY")
PREDICTED PRICE =
> mpredict("SPY")
PREDICTED PRICE =
> mpredict("SPY")
PREDICTED PRICE =
> mpredict("SPY")
PREDICTED PRICE =
135.4648 ACTUAL PRICE = 132.86
124.4090 ACTUAL PRICE = 132.86
133.5763 ACTUAL PRICE = 132.86
128.8359 ACTUAL PRICE = 132.86
136.2620 ACTUAL PRICE = 132.86
52
Dont even think about AAPL
Ugh!
> mpredict("AAPL")
PREDICTED PRICE =
> mpredict("AAPL")
PREDICTED PRICE =
> mpredict("AAPL")
PREDICTED PRICE =
> mpredict("AAPL")
PREDICTED PRICE =
> mpredict("AAPL")
PREDICTED PRICE =
390.8149 ACTUAL PRICE = 335.06
359.5131 ACTUAL PRICE = 335.06
358.1080 ACTUAL PRICE = 335.06
434.0794 ACTUAL PRICE = 335.06
330.2102 ACTUAL PRICE = 335.06
53
What if we average lots of predictions?
Still basically sucks
> mpredictavg("GE")
AVG PREDICTED PRICE =
> mpredictavg("AAPL")
AVG PREDICTED PRICE =
> mpredictavg("SPY")
AVG PREDICTED PRICE =
> mpredictavg("SPY")
AVG PREDICTED PRICE =
> mpredictavg("JNJ")
AVG PREDICTED PRICE =
> mpredictavg("SSO")
AVG PREDICTED PRICE =
> mpredictavg("SDS")
AVG PREDICTED PRICE =
> mpredictavg("DIA")
AVG PREDICTED PRICE =
>
20.89974 ACTUAL PRICE = 20.19
363.9570 ACTUAL PRICE = 335.06
132.3680 ACTUAL PRICE = 132.86
132.4612 ACTUAL PRICE = 132.86
59.39538 ACTUAL PRICE = 59.46
52.73392 ACTUAL PRICE = 53.47
21.18263 ACTUAL PRICE = 20.83
121.6629 ACTUAL PRICE = 123.65
Keep in mind, were focusing on predicting one day
54
Hmmmm
Predicting prices is hard
What if all we want to do is predict direction?
Hmmm
That is, at the start the 20 day period, we
want to know whether to go long or short the
stock.
Would our prediction work as a directional
indicator?
55
More R code (phew)
Code fragment
if (as.numeric(last(last20))>as.numeric(last(Cl(sdata)))) tdir="UP"
if (as.numeric(last(last20))<as.numeric(last(Cl(sdata)))) tdir="DOWN"
if (as.numeric(last(last20))==as.numeric(last(Cl(sdata)))) tdir="FLAT"
if (as.numeric(mean(snewp))>as.numeric(last(Cl(sdata)))) pdir="UP"
if (as.numeric(mean(snewp))<as.numeric(last(Cl(sdata)))) pdir="DOWN"
if (as.numeric(mean(snewp))==as.numeric(last(Cl(sdata)))) pdir="FLAT"
56
Hmm-does it work or not?
> mdirection("GE")
----------- STOCK = GE ----------------------Day 21 = 20.36 PREDICTION = 20.86483 LAST DAY = 20.19
PREDICTED DIR = UP ACTUAL DIR = DOWN
--------------------------------> mdirection("AAPL")
----------- STOCK = AAPL ----------------------Day 21 = 351.99 PREDICTION = 364.0584 LAST DAY = 335.06
PREDICTED DIR = UP ACTUAL DIR = DOWN
--------------------------------> mdirection("SPY")
----------- STOCK = SPY ----------------------Day 21 = 130.84 PREDICTION = 132.3852 LAST DAY = 132.86
PREDICTED DIR = UP ACTUAL DIR = UP
--------------------------------> mdirection("DIA")
----------- STOCK = DIA ----------------------Day 21 = 120.42 PREDICTION = 121.5847 LAST DAY = 123.65
PREDICTED DIR = UP ACTUAL DIR = UP
--------------------------------> mdirection("SDS")
----------- STOCK = SDS ----------------------Day 21 = 21.75 PREDICTION = 21.09891 LAST DAY = 20.83
PREDICTED DIR = DOWN ACTUAL DIR = DOWN
---------------------------------
57
Try the NASDAQ 100 stocks
----------- STOCK = ALTR ----------------------Day 21 = 40.46 PREDICTION = 42.10049 LAST DAY = 42.81
PREDICTED DIR = UP ACTUAL DIR = UP
------------------------------------------- STOCK = AMZN ----------------------Day 21 = 168.07 PREDICTION = 170.7999 LAST DAY = 184.71
PREDICTED DIR = UP ACTUAL DIR = UP
------------------------------------------- STOCK = AMGN ----------------------Day 21 = 53.53 PREDICTION = 53.54674 LAST DAY = 53.9
PREDICTED DIR = UP ACTUAL DIR = UP
------------------------------------------- STOCK = APOL ----------------------Day 21 = 42.29 PREDICTION = 41.51876 LAST DAY = 42.11
PREDICTED DIR = DOWN ACTUAL DIR = DOWN
------------------------------------------- STOCK = AAPL ----------------------Day 21 = 351.99 PREDICTION = 364.7487 LAST DAY = 335.06
PREDICTED DIR = UP ACTUAL DIR = DOWN
---------------------------------
58
Overall Statistics
When we apply this technique to the nasdaq
100 stocks (for one particular day! Not very
scientific).
It is accurate on 67% of the stocks.
Better than flipping a coin.
To really test this we would have to try
numerous days (and vary the number of days
to predict for, and to use for the density
estimate and this is why hedge funds hire
statisticians.)
59