But a pleasant 3 days of sleep deprivation it was — seeing old friends and making new ones.

I’m not sure that Mother Mary believed me that in our house we still use the mug from the original R in Finance.  Here’s proof that it still survives when so many of its comrades have perished.

I had great fun giving my presentation.  The annotated slides and source for the R package are available on the Burns Statistics website.

Highlights

There were lots of good parts.  I’ll limit myself to three (three is the largest number of chunks that people can easily remember).

First

Ahead by 1609.344 meters is Doug Service.  Universally my experience with talks involving stochastic calculus is that my eyes glaze over before the first dW drops.  At the start of Doug’s presentation my eyes were doing their cross-check preparing for landing.  But the landing never happened.  I have no real idea of what the equations were about — that didn’t matter: Doug made the idea of what he was doing perfectly clear.

If you insist on having formulas in your presentation (I highly discourage you from such a point of view), then you must track Doug down and make him teach you how to do it right.

Second

Second is Dirk’s update on the package interfacing Bloomberg with R:

Last year I said A.

So I lied.

Actually …

Last year I said B.

So I lied.

Actually …

Last year I said C.

So I lied.

Actually …

(You should be starting to see a pattern.)

Third

Tied for third are:

Matt Dziubinski’s high performance talk.  High performance was the topic and high performance was the delivery.  Clearly Matt has at least four cores dedicated to his mouth, with negligible latency.

Bryan Lewis’ politics of big data hardware implementation structures.

In the end

The taxi back to the airport was flanked by a full moon on the left and the rising sun on the right.

I was born to rock the boat
Some may sink but we will float
— from “Mutineer” by Warren Zevon

R in Finance

2016 May 20-21, Chicago.

2 days, limited space, 50 speakers, including:

Pat Burns on “Some Linguistics of Quantitative Finance”

Abstract: How can the abstract be written for a talk with an ambiguous and possibly misleading title without itself being vague and misleading? I don’t know, but perhaps: A quest to discover how markets work that starts in the Yucatan.

See the conference website.

Other Events

London Quant Group Spring Seminar

2016 May 12, London.

See the LQG website for details.

Forecasting Financial Markets

2016 May 25-27. Hanover.

See the conference website for details.

 UseR! 2016

2016 June 27-30. Stanford.  Apparently sold out.

The UseR! 2016 website.

London Quant Group Autumn Seminar

2016 September 11-14, Oxford. 30th annual celebration event.

See the LQG website for details.

EARL 2016

2016 September 13-15, London. Effective Applications of the R Language.
See the EARL conference website for details.

Even more events

MoneyScience has an events calendar.

Executive summary

If you are confused by economics, then you should read this book.  If you are not confused by economics, then you really, really need to read this book.

The book is inexpensive, not very long, and it’s important — there is no reason not to read it and good reason to read it (including some splashes of humor).

Structure

The book starts by explaining four seemingly random scientific revolutions:

• The earth goes around the sun, rather than vice versa
• Blood circulates in the body
• Natural selection via survival of the fittest
• Continental drift and plate tectonics

By the end of the book we see that this isn’t such a random selection of scientific issues after all.

Then there is a chapter that explains the viewpoints of the various economic schools.  This is the clearest explanation of the economic camps that I’ve ever seen — it erased several points of confusion for me.

Finally, a synthesis of the warring factions is proposed that explains many of the anomalies.  George calls this the circulatory growth model.  It must be a good model — there’s something in it for pretty much everyone to hate.

Blindness

People tend to ignore what their theories don’t explain.  Economists have learned this well:

• luxury goods should not exist
• I should always buy the cheapest of anything
• the size of my neighbor’s television shouldn’t affect the size of television I want
• economic growth should be more randomly scattered through time, and location

But there is one thing that is ignored in spectacular fashion — government spending. Many economists act as if government doesn’t exist at all.  There’s a disease like that:

Patients with unilateral neglect behave not only as if nothing were actually happening in the left hemispace, but also as if nothing of any importance could be expected to occur there.

— from The Man Who Mistook His Wife For A Hat by Oliver Sacks

Quotes

It appears that we humans cannot accept having no theory … we cannot be persuaded to step from a wrong theory into the abyss while we wait for a better theory to come along.  This is why, regardless of how many times the behaviour of the economy refutes current economic theory, that theory will not be rejected until a better model is proposed.

… small government, light regulation and low taxation is the economic condition that has persisted across most of the world through all of history.  Why then has the world not been plagued by a permanent growth-boom lasting thousands of years?

In 2008… everyone suddenly rejected the silly notion of an inherently stable economy and embraced Minsky’s idea of an inherently unstable economy … In the years since 2008, Minsky has been forgotten and the herd has quickly drifted back to recolonise the neoclassical zone, conveniently forgetting that it ever left.

… to understand and describe [our current economy] cannot be achieved by building abstract models of imaginary alternate economies, where governments are absent, markets are perfect and people act as machines rather than humans.

Favorite sentence

With democratic capitalism we have inadvertently stumbled upon a system that lets us exploit the best of our competitive nature while managing the worst of our selfish nature.  We should celebrate this triumph.

Quibbles

I think that the dissipation of the happiness of lottery winners is a more complex phenomenon than the book implies.

Great as it is, I think there is something missing from the circulatory growth model.  I’m not sure what is missing though.

Related books

I’ve reviewed The Origin of Financial Crises, George’s previous book.  The books are related, but I don’t see that it matters which you read first.

A great book that is surprisingly not referenced is The Truth about Markets by John Kay.  This supports the argument that suitable institutions are necessary for economic growth.

Follow-up action

Kuhn’s final blow to the ego of scientists was his observation that the vanguard of the new paradigm is usually led by non-expert outsiders.  These outsiders are often thoroughly frustrated with the confusion of the old order and eager to accept the clarity of a new idea.

If Kuhn is correct, then we can’t leave it to economists to sort themselves out.  The stakes of getting economics right can’t be higher.

Suggested next steps:

• Read George’s blog
• Follow George on Twitter: @george_cooper__
• Tell your friends

Epilogue

Brian: You’ve got to think for yourselves.  You’re all individuals.

Crowd (in unison): Yes, we’re all individuals.

Brian: You’re all different.

Crowd (in unison): Yes, we’re all different.

Lone voice: I’m not.

from The Life of Brian by Monty Python

The first EARL Conference (Effective Applications of the R Language) was held 2014 September 15-17 in London.

Talk

My talk was “Effective risk management with R” (annotated slides).

Instability hypothesis

When I was preparing for the talk, one of my ideas was to show the Google trend for searches for Minsky’s instability hypothesis.  I thought it would look a lot like Figure 1.

Figure 1: My guess of interest in the instability hypothesis over time. My impression is that the instability hypothesis is pretty much totally ignored except immediately following market crashes.  That’s the sort of thinking that allows the pattern to continue.

Figure 2: Google trend for: minsky instability hypothesis.

The trend that I actually saw was severely different from my prior.  This leads me to pretty much entirely discount the data.  Part of the problem is the small number of searches at any time.

Epilogue

I’ll tell ya truly that I sometimes lie

— from “What Kinda Guy?” by Steve Forbert

EARL

As in “Effective Applications of the R Language”.

2014 September 15-17, London.

Somehow they gave higher billing to Ben Goldacre than to Pat Burns.  If Obama were coming, they’d probably bill him above me too — and what does he know about R?  In spite of that little glitch, I’m sure it will be a good time.

See the conference website.

New Events

QWAFAFEW — New York

2014 August 26.  Correcting quantitative finance’s erroneous understanding and calculations of mathematical averaging

London Quant Group Autumn Seminar

2014 September 7-10, Oxford.

See the LQG website for details.

14-10 Club

2014 September 11.

“Bitcoin: a human reaction to centralisation” and “Happiness by Design”.

Details at the 14-10 website.

London Quant Group

2014 October 7, London.

See the LQG website for details.

14-10 Club

2014 October 9.

“Advanced internal funding process for banks” and “Behavioural Finance”.

Details at the 14-10 website.

14-10 Club

2014 November 6.

“What can science teach business (and vice versa)” and “Physics of Finance”.

Details at the 14-10 website.

London Quant Group

2014 November 11, London.

See the LQG website for details.

London Quant Group

2014 December 8, London.

See the LQG website for details.

14-10 Club

2014 December 10.

“Hedge funds turn to psychology software” and “The softer side of physics”.

Details at the 14-10 website.

Computational and Financial Econometrics 2014

2014 December 6 – 8, Pisa.

Details at the conference website.

Even more events

MoneyScience has an events calendar.

Variance estimation

The major functionality in the package is variance estimation:

• Ledoit-Wolf shrinkage via var.shrink.eqcor
• statistical factor model (principal components) via factor.model.stat

There have been a number of previous blog posts on both factor models and Ledoit-Wolf shrinkage.

Positive-definiteness

The default value of the tol argument to var.shrink.eqcor is different in the new version.  The blog post “Correlations and positive-definiteness” suggested that the previous default was too extreme.

In finance we almost always want a variance matrix to definitely be positive-definite.  This is especially true when doing optimization.  Positive-definite means that all eigenvalues are greater than zero.  The tol argument in var.shrink.eqcor gives the minimum size of an eigenvalue in the result relative to the largest eigenvalue.

An eigenvalue of zero means that the variance matrix is saying that there is a portfolio (described by the corresponding eigenvector) that has zero variance.  An optimizer is going to really, really like such a portfolio — who wouldn’t?

The previous default for tol was 0.001.  The new default is 0.0001.  This is still a guess, but hopefully a more appropriate guess.

Efficiency

Both the input matrix and the resulting variance matrix can be large.  Financial applications often have a few thousand assets in a variance.  If there are multiple copies of these large objects within the function evaluation, that can be enough to push memory use over the limit in a computation.

Originally var.shrink.eqcor and factor.model.stat were written to get the computations right rather than worrying about memory or time efficiency.  (That’s the right order in which to worry.)  The new version uses fewer copies of the large objects.  That change accidentally sped up factor.model.stat a little.

If you run out of memory using one of these two functions, then:

First consider moving to 64-bit if you are still using 32-bit — R runs the same on either.  If you are on 64-bit already, consider getting more RAM (an application of Uwe’s maxim).

Next, you could try to make your code more efficient.  In particular, avoid Circle 2 of The R Inferno.

If none of that is applicable, then you may want to take evasive action.  var.shrink.eqcor uses more memory than factor.model.stat.  If you run out of memory with factor.model.stat, then you can set iter.max=0 to save some memory (and time) at the expense of a (probably) slightly worse estimate.

 Weird data

Real data (financial and other) tend not to be nice textbook examples.  Often there are missing values.  In finance we even want to have variables (assets) in a variance matrix that have no data at all.

One of many possible data glitches in finance is for a price to get stuck on one value.  This translates to returns that are all zero.  Zero is an underestimate of the variance of an asset, and most likely suggests there is a data problem.  But that isn’t true of all settings.  factor.model.stat has an argument called constant.returns.okay.  The new version produces a different — that is, not backward-compatible — result when this argument is TRUE (which is not the default).

It is reasonable to view the previous behavior regarding this argument to be a bug:  All correlations involving the variable are zero and the variance is a certain quantile of the variances of the other variables.  The new behavior is to set the variance to zero as well.

Probably the best thing to do when confronted with this is to set the constant returns to missing values.

A possible self-inflicted weirdness is to give the variance estimators prices rather than returns.  By default the functions produce a warning if the input data are all non-negative.

Time weights

Financial markets are dynamic.  That means that old data goes stale.  The most typical approach to this problem is to equally weight some number of the most recent observations.  The default weighting in var.shrink.eqcor and factor.model.stat is to linearly decrease weights in such a way that the oldest observation given has one-third the weight of the most recent.

The new function in the package, called slideWeight, allows both of these, plus compromises between them.  The typical shape of weights from this function is shown in Figure 1.

Figure 1: Slide weights where 50 observations have full weight and 170 have some weight. 

The slideWeight function lets you specify a certain number of observations to have full weight, and a value beyond which all weights are zero.  The weights fall linearly between these two points.

See also

If you are new to R, learning R may be easier than you think.

Simple question

If I have correlation matrices each estimated with a month of daily returns, how much worse is the average of six of those compared to the estimate with six months of daily data?

Expected answer

Do a statistical bootstrap with the returns and compare the standard deviations across bootstrap samples for each correlation.  This will show how much bigger the monthly standard deviations are compared to the daily standard deviations.

Return data

We use 126 daily (essentially 2012 H2) returns of 100 large cap US stocks.  Figure 1 is a scatter plot of the 4950 correlations estimated with the two methods.

Figure 1: Correlations of returns estimated by mean monthly versus daily. 

Surprise #1: Figure 1 is not boring.

Figure 2 compares the standard deviations of bootstraps of the two estimation methods.

Figure 2: Standard deviations across bootstrap samples of the monthly estimation versus daily estimation. 

Surprise #2: The wrong estimator tends to be more variable.

Random normal data

We can experiment:  Create new data shaped like our original data but that is randomly generated as multivariate normal with correlation equal to the daily estimate.

Figure 3 is like Figure 1 but for the normal data.

Figure 3: Correlation estimates: monthly versus daily for the randomly generated normal data. 

Semi-surprise #1: Figure 3 is boring, unlike Figure 1.

We know that the bootstrap is appropriate for the normal data: the observations are independent and identically distributed — perfect for the bootstrap.  Figure 4 is the normal data version of Figure 2.

Figure 4: Standard deviations across bootstrap samples of the monthly estimation versus daily estimation of the randomly generated normal data. 

Now that we have fake data for which we actually know the true value of the correlations, we can look at the mean squared error of the estimates.  This is done in Figure 5.

Figure 5: Mean squared errors (MSE) across bootstrap samples of the monthly estimation versus daily estimation of the randomly generated normal data. 

Surprise #3: There’s hardly any difference at all between the quality of the two estimates for the normal data.

Back to reality

We can highlight the correlations where the two methods differ significantly, as shown in Figure 6.

Figure 6: Correlations estimated by mean monthly versus daily with strange points highlighted (in magenta). 

Figure 7 shows where those highlighted points show up in the plot of bootstrap standard deviations.

Figure 7: Standard deviations across bootstrap samples of the monthly estimation versus daily estimation with correlation outliers highlighted. 

Surprise #4: The disagreement over negative daily correlations does not explain the relatively low variability of the monthly estimates — at least not entirely.

Note that both standard deviations for the real data tend to be larger than for the normal data.

We can define weird the other way around: define extreme in terms of the bootstrap standard deviations and then see where those show up in the correlation scatterplot.  This is Figures 8 and 9.

Figure 8: Standard deviations across bootstrap samples of the monthly estimation versus daily estimation with standard deviation outliers highlighted. 

Figure 9: Correlations estimated by mean monthly versus daily with standard deviation outliers highlighted. 

Summary

I don’t know what is going on.  There is some sort of dependence that is showing up in the correlations of the returns.  Part of the puzzle is probably that return correlations (and means) are not constant through time.

Epilogue

Complications in the air
Complications in the forecast now
Complication everywhere

from “Complications” by Steve Forbert

Appendix R

The computations were in R.

daily correlation estimate

Getting the correlation matrix using daily data is easy enough:

corDay <- cor(retmat6m)

mean monthly estimate

A little more goes into our monthly estimate.  The first thing we do is create a three dimensional array to hold the six monthly estimates:

corMons <- array(NA, c(100, 100, 6), 
   c(dimnames(corDay), list(NULL)))

Note that c is used to put two lists together.

Then we compute the correlations in each month:

dseq <- 1:21
for(i in 1:6) {
   corMons[,,i] <- cor(retmat6m[dseq + (i-1) * 21,])
}

Finally, we average the correlations, ending up with a matrix:

corMM <- apply(corMons, 1:2, mean)

(There is a chapter in S Poetry about higher dimensional arrays that may be of use if you start working with them.)

scatter plot

We just want one copy of the off-diagonal values in the correlation matrix.  We can get a logical vector selecting the lower triangle with:

clt <- lower.tri(corDay, diag=FALSE)

(The default value of diag is FALSE, but I typically include diag in my calls to lower.tri just to make sure I have it the right way around without looking.)

This is used like:

plot(corDay[clt], corMM[clt])

write a simple function

We can write a little function that does the mean monthly correlation estimate:

pp.cor6m <- function(x)
{
   mc <- array(NA, c(ncol(x), ncol(x), 6))
   dseq <- 1:21
   for(i in 1:6) {
      mc[,,i] <- cor(x[dseq + (i-1) * 21,])
   }
   apply(mc, 1:2, mean)
}

bootstrap correlations

We start off again by creating arrays that we will then populate:

bootcord <- bootcorm <- array(NA, c(100, 100, 200))

Now do the actual bootstrapping:

for(i in 1:200) {
   bx <- retmat6m[sample(126, 126, replace=TRUE),]
   bootcord[,,i] <- cor(bx)
   bootcorm[,,i] <- pp.cor6m(bx)
}

Finally, get the standard deviation across the bootstrap estimates for each correlation:

bootd.sd <- apply(bootcord, 1:2, sd)
bootm.sd <- apply(bootcorm, 1:2, sd)

generate normal data

There is a function in the MASS package for generating random multivariate normal data.

require(MASS)

This is used like:

normmat <- mvrnorm(126, mu=rep(0,100), Sigma=corDay)

bootstrap MSE

The bootstrapping is:

nbtcord <- nbtcorm <- array(NA, c(100, 100, 200))
for(i in 1:200) {
   bx <- normmat[sample(126, 126, replace=TRUE),]
   nbtcord[,,i] <- cor(bx)
   nbtcorm[,,i] <- pp.cor6m(bx)
}

Then the mean squared errors are calculated like:

nbtd.mse <- apply(nbtcord - as.vector(corDay), 1:2, 
                  function(x) mean(x^2))

The as.vector is used because there is an error otherwise. When arithmetic is done on arrays in R, there is a check to see if they are conformable. A matrix and a three dimensional array are not conformable.  Removing the matrixness of corDay allows the computation to proceed, and it is the right computation in this case.

define strange

When looking at Figure 1 we have a sense of which points are weird.  There are multiple ways to translate our intuition into a definition.  Here is the one I used:

corMM[clt] - corDay[clt] > .1 & corDay[clt] < 0

That is, the monthly estimate needs to be at least 0.1 bigger than the daily estimate and the daily estimate must be negative.

The definitions of weird in the other direction are:

wei <- bootd.sd[clt] - bootm.sd[clt] > .02
wei2 <- bootd.sd[clt] - bootm.sd[clt] > .04

Previously

Problem 4 of “The top 7 portfolio optimization problems” concerns the use of higher moments.

“Further adventures with higher moments” is the most recent in a series of posts on the efficacy of higher moments in optimization.  This set includes the observation that “trade selection” is a better term than “portfolio optimization”.

Choices

For any particular optimization problem there are a few possibilities of what you might believe about higher moments:

1. The returns are normally distributed
2. Adding higher moments to the optimization has no benefit
3. The effect of adding higher moments is unknown
4. Can’t be bothered to add higher moments even though they would help
5. Use higher moments

1) Normally distributed returns

If you believe this, then you might miss the satire in “The distribution of financial returns made simple”.

In this case higher moments are known, and of no use.

2) No benefit of higher moments

There is persistent misinformation that this possibility doesn’t exist.  The commonly repeated message is that mean-variance optimization implies that you believe returns are normal.

That is not the case.  It has things backwards: if returns are normal, then any optimization is mean-variance (though it could be done inefficiently by adding higher moments).

Doing mean-variance optimization is perfectly okay when the returns have symmetric distributions.  It is the best thing to do if there is no predictability of higher moments.

3) Unknown effect

This is appropriate, for instance, when it is clear that there is some predictability of skewness and kurtosis, but it is not clear that the predictability will translate into better results.

4) Can’t be bothered

Incorporating higher moments would improve the utility of the portfolio, but seemingly not the overall utility that includes the process of managing the portfolio.

5) Include higher moments

Higher moments are put into the optimization — either explicitly, or implicitly via a scenario optimization.

Credo

For most assets other than equities, my bias would be to do scenario optimization.  High-level asset allocation would generally be an exception.

For equities I am somewhere among beliefs 2 through 4.  Before I started looking at the stability of skewness and kurtosis among random portfolios, I was leaning towards 4.  Since I have done that experimenting, I am much closer to 2.

Epilogue

I believe the stars keep shining all through the night.
I believe if we just keep trying it will be all right.

from “I Believe” by Chris Isaak

Definition

Volatility is the annualized standard deviation of returns — it is often expressed in percent.

A volatility of 20 means that there is about a one-third probability that an asset’s price a year from now will have fallen or risen by more than 20% from its present value.

In R the computation, given a series of daily prices, looks like:

sqrt(252) * sd(diff(log(priceSeriesDaily))) * 100

Usually — as here — log returns are used (though it is unlikely to make much difference).

Historical estimation

What frequency of returns should be used when estimating volatility?  The main choices are daily, weekly and monthly.  (In some circumstances the choices might be something like 1-minute, 5-minute, 10-minute, 15-minute.)

There is folklore that it is better to use monthly data than daily data because daily data is more noisy.  But noise is what we want to measure — daily data produces less variable estimates of volatility.

However, this is finance, so things aren’t that easy.  There is the phenomenon of variance compression — that is, it isn’t clear that volatility from monthly and daily data are aiming at the same target.

Another complication is if there are assets from around the globe.  In that case asynchrony is a problem with daily data.  For example, the daily return of assets traded in Tokyo, Frankfurt and New York for January 21 are returns over three different time periods.  That matters (not much for volatility, but definitely for correlations).  Weekly returns as well as fancy estimation methods alleviate the problem.

Through time

Volatility would be more boring if finance were like other fields where standard deviations never change.

Financial markets seem to universally experience volatility clustering.  Meaning that there are time periods of high volatility and periods of (relatively) low volatility.

But why?

Changes in volatility would make sense if there were periods of greater uncertainty about the future.  But the future is always uncertain.

The day after Lehman Brothers failed there was a great deal more uncertainty about the markets than the week before.  However, that uncertainty was generated internally to the markets, it wasn’t greater uncertainty about the physical world.

Agent-based models of markets easily reproduce volatility clustering.  Rather than agents changing their minds about the best strategy at random times, mind-changing tends to cluster.  When lots of agents change strategies, then the market volatility goes up.

Across assets

Different assets in the same market have different volatilities.  But why?

Some assets are especially responsive to the behavior of the  rest of the market.  Some assets are inherently more uncertain.

Factor models explicitly divide these two ideas into systematic and idiosyncratic risk.

Implied volatility

Implied estimates come from using a model backwards.  Suppose you have a model that approximates X by using variables A, B and C.  If you know X, A and B, then you can get an implied estimate of C by plugging X, A and B into the model.

An example is implied alpha.

Implied volatility comes from options pricing models.  The most famous of these is the Black-Scholes formula.

Black-Sholes uses some assumptions and derives a formula to say what the price of an option should be.  The items in the formula are all known except for the average volatility between the present and the expiration date of the option.

We can use known options prices and the Black-Scholes formula to figure out what volatility is implied.  You can think of this as the market prediction of volatility — assuming the model is correct (but it isn’t).

Not risk

Some people equate volatility with risk.  Volatility is not risk.

Volatility is one measure of risk.  In some situations it is a pretty good measure of risk.  But always with risk there are some things that we can measure and some things that escape our measurement.

Part of risk measurement should be trying to assess what is left outside.

Epilogue

Before I built a wall I’d ask to know 
What I was walling in or walling out, 
And to whom I was like to give offense.

from “Mending Wall” by Robert Frost

Appendix R

I made the statement that using simple returns rather than log returns is unlikely to make much difference.  We shouldn’t take my word for it.  We can experiment.

In most blog posts I give only a sample of the R code.  Here, though, is the complete R code for this experiment.

write a function

pp.dualvol <- function(x, frequency=252) {
   lv <- sd(x)
   sv <- sd(exp(x)) # same as sd(exp(x) - 1)
   100 * sqrt(frequency) * 
      c(logReturn=lv, simpleReturn=sv)
}

get starting points of the samples

starts <- sample(length(spxret) - 253, 1000)

This statement presumes that the number sampled is smaller than the number of values being sampled since it is using the default of sampling without replacement.

compute volatilities

volpairs <- array(NA, c(1000, 2), 
   list(NULL, names(pp.dualvol(tail(spxret, 10)))))

tseq <- 0:62

for(i in seq_along(starts)) {
   volpairs[i,] <- pp.dualvol(spxret[starts[i] + tseq])
}

write plot function

P.dailyVolPair63 <-
 function (filename = "dailyVolPair63.png") 
 {
    if(length(filename)) {
       png(file=filename, width=512)
       par(mar=c(5,4, 0, 2) + .1)
    }
    plot(volpairs, col="steelblue")
    abline(0, 1, col="gold")
    if(length(filename)) {
       dev.off()
    }
 }

do the plot

P.dailyVolPair63()

Figure 1: Volatility estimates on a quarter of daily data — simple versus log returns. 

We can see that almost all of the 1000 estimates are virtually the same when using either type of return.  There are a few large volatilities that are slightly different.  The differences will be due to large returns, which means that returns at lower frequencies will yield bigger differences between the volatility estimates.

Previously

There have been several posts about garch.  In particular:

• A practical introduction to garch modeling
• The components garch model in the rugarch package

Genesis

A reader emailed me because he was confused about the workings of garch in general, and simulation with the empirical distribution in particular.

If there is one confused person, there are more.

Model view

A basic feature of garch models is that they are in discrete time.  The frequency is usually daily, in which case we pretend that the volatility of the returns is constant throughout each day.

Deep in the heart of a garch model is an innovation at each timepoint.  This is — conceptually — a draw from a statistical distribution that has variance 1.  Note that variance 1 rules out infinite variance, and hence stable distributions.  Also, the normal distribution is unlikely to be a good approximation.

The innovation is the net result of the “news” that comes into the market during the period.  If the innovation is large (in absolute value, and compared to a variance 1 distribution), then that is a signal that future volatility is likely to be larger.

But we don’t see innovations.  We see innovations scaled by the current volatility.

What the garch model keeps track of is the conditional variance at each time.  The “conditional” being short for conditional on the returns (and possibly other variables) that have already occurred.  “Variance” because it is mathematically handy — in higher dimensions there are variance matrices but you’ve never heard of a standard deviation matrix.

So far our model is:

return at time t  =  innovation at time t, scaled using the conditional variance at time t

There is one more element to the model: the mean return at time t.  The typical garch model is:

return at time t  =  mean return at time t  +  innovation at time t, scaled using the conditional variance at time t

Keep in mind that the conditional variance needs to be transformed before it can be used to scale the innovation.

Reality view

The fact is that there is a return — a single number — for a time period.  The garch models are taking that number and decomposing it into three pieces:

• the mean return at that time
• the conditional variance at that time
• the innovation for that time

Those three things are fictional.  And the fictions of the last two heavily depend on each other.

Sometimes they are useful fictions.

Data view

The basic input to a garch model is a time series of returns.  It is log returns that make the most sense in this setting.

Possible specifications for the mean return are:

• zero
• a constant (to be estimated)
• a simple autoregression — AR(1)
• an ARMA model
• something even fancier

There is unlikely to be much difference between assuming the mean is zero and any of the others.  In particular, the estimation of the variance is highly unlikely to change significantly.

Residuals

It is easy for the term “residual” to be confusing in garch models.  That’s because the residuals don’t really have anything to do with the garchiness of the model.

The residual at time t is the return at time t minus the estimated mean return at time t.

More garchy — and more useful — are the standardized residuals.  The standardized residual at time t is the residual at time t divided by the square root of the conditional variance at time t.

In other words, the standardized residuals are the estimates of the innovations.

Fit

There are two parts to the fit of a garch model:

• fit for the mean return
• fit for the conditional variance (or standard deviation or volatility)

The second of these is the more interesting.

Simulation

To do a garch simulation you need to start with:

• a garch model with the parameters specified
• a state for the conditional variance

Usually, of course, the parameter specification is taken from the estimated model.  It is also typical that the variance state is the most recent variance state from the estimated model.

There is one more thing the simulation needs — innovations.

The two most common choices for generating a series of innovations are:

• sample from the distribution assumed in the model (including any estimated parameters)
• sample from the standardized residuals of the estimated model

The process of simulation is to use the first innovation, scale it via the volatility state and add in the mean return to get the first return in the simulation — let’s say this is for time T+1.

We then need to update the volatility state for T+1 — this depends on the innovation.

After those two steps we are in the same situation except we are one time step further.  We continue doing those steps until we have as many timepoints as we want.

Usually we do that process multiple times so that we have a large number of simulated series.

Epilogue

Your brow is sweatin’ and your mouth gets dry,
Fancy people go driftin’ by.
The moment of truth is right at hand,
Just one more nightmare you can stand.

from “Stage Fright” by Robbie Robertson

Appendix R

You can use garch in R, among other places.

There are alternatives, but the rugarch package seems to be the most complete implementation in R.  Here is a small example.  The first step is:

require(rugarch)

This command will only work once the package has been installed on your machine.

set the garch specification

We start by creating a specification of the garch model:

comtspec <- ugarchspec(mean.model=list(
  armaOrder=c(0,0)), distribution="std",
  variance.model=list(model="csGARCH"))

This specification does three things. It sets:

• the model for the mean return to be a constant to be estimated — the default is an ARMA(1,1)
• the distribution of the innovations to be a Student-t distribution — the default is normal
• the model for the variance to be components garch — the default is garch(1,1)

estimate the model

We use the specification when fitting the model:

garEstim <- ugarchfit(comtspec, data=tail(spxret, 2000))

The data are the most recent 2000 returns in the vector.

get residuals

The regular residuals for the garch model are:

garResid <- residuals(garEstim)

The standardized residuals are computed like:

garStandResid <- residuals(garEstim, standardize=TRUE)

We can see that the results make sense by inspecting the variances:

> var(tail(spxret, 2000))
[1] 0.0001969673
> var(garResid)
             [,1]
[1,] 0.0001969673
> var(garStandResid)
         [,1]
[1,] 1.204764

I’m slightly concerned that the variance of the standardized residuals is so far from 1, but clearly it is trying to be close to 1.

how you might have found how to get residuals

The residuals function is generic.  There are two ways to look for the methods of a generic function:

> methods('residuals')
[1] residuals.default        residuals.glm 
[3] residuals.HoltWinters*   residuals.isoreg* 
[5] residuals.lm             residuals.nls* 
[7] residuals.smooth.spline* residuals.tukeyline*
   Non-visible functions are asterisked

> showMethods('residuals')
Function: residuals (package stats)
object="ANY"
object="ARFIMAfilter"
object="ARFIMAfit"
object="ARFIMAmultifilter"
object="ARFIMAmultifit"
object="uGARCHfilter"
object="uGARCHfit"
object="uGARCHmultifilter"
object="uGARCHmultifit"

(Note that the results you get depend on what packages are loaded into your R session.)

The class of the object that we want is:

> class(garEstim)
[1] "uGARCHfit"
attr(,"package")
[1] "rugarch"

This is found in the second group.  To have a look at the method in this group, you do:

getMethod('residuals', 'uGARCHfit')

It is here that we see that standardize is an argument (though not exactly blatantly obviously).

encore

The post “garch and the distribution of returns” has some examples of simulation.