The following content is provided
under a Creative Commons license.
Your support will help MIT OpenCourseWare continue
to offer high quality educational resources for free.
To make a donation or view additional materials
from hundreds of MIT courses, visit MIT OpenCourseWare
at ocw.mit.edu.
PROFESSOR: All right.
I want to complete the discussion
on volatility modeling in the first part of the lecture
today.
And last time we addressed the definition
of ARCH models, which allow for time varying
volatility in modeling the returns of a financial time
series.
And we were looking last time at modeling the Euro-dollar
exchange rate returns.
And we went through fitting arch models to those returns,
and also looked at fitting the GARCH model to those returns.
And to recap, the GARCH model extends upon the ARCH model
by adding some extra terms.
So if you look at this expression for the GARCH model,
the first two terms for the time varying volatility sigma
squared t is a linear combination
of the past sort of residual returns squared.
That's the ARCH model, p [? of ?] [? those. ?] So
the current volatility depends upon what's happened in excess
returns over the last p periods.
But then we add an extra term, which
is corresponds to queue levels of the previous volatility.
And so what we're doing with GARCH
models is adding extra parameters to the ARCH,
but an advantage of considering these extra parameters which
relate basically the current volatility sigma squared t
with the previous or lagged value sigma squared t
minus j for lags j is that we may
be able to have a model with many fewer parameters.
So indeed, if we fit these models to the exchange rate
returns, what we found last time-- let me go through
and show that-- was basically here
are various fits of the three cases of ARCH models.
ARCH orders 1, 2, and 10, thinking we maybe
need many lags to fit volatility.
And then the GARCH model 1,1, where we only
have one ARCH term and one GARCH term.
And so basically the green line-- or rather the blue line
in this graph, shows the plot of the fitted GARCH 1,1,
model as compared with the ARCH models.
Now, in looking at this graph, one
can actually see some features of how these models are
fitting volatility, which is important to understand.
One is that the ARCH models have a hard lower
bound on the volatility.
Basically there's a constant term
in the volatility equation.
And because the additional terms are squared, excess returns,
it basically-- the volatility does
have lower bound of that intercept.
So depending on what range you fit the data over,
that lower bound is going to be defined by--
or it will be determined by the data you're fitting to.
As you increase the ARCH order, you basically
allow for a greater range of-- or a lower lower bound of that.
And with the GARCH model you can see
that this blue line is actually predicting
very different levels of volatility
over the entire range of the series.
So it really is much more flexible.
Now-- and in these fits, we are assuming Gaussian distributions
for the innovations in the return series.
We'll soon pursue looking at alternatives to that,
but let me talk just a little bit more
about the GARCH model going back to lecture notes here.
So let me expand this.
OK.
So there's the specification.
The GARCH 1,1 model.
One thing to note is that this GARCH 1,1 model does
relate to an ARMA, an Auto Regressive Moving Average,
process in the squared residuals.
So if we look at the top line, which
is the equation for the GARCH 1,1 model,
consider eliminating sigma squared t
by using a new innovation term, little ut, which
is the difference between the squared residual
and the true volatility given by the model.
So if you plug in the difference between our squared excess
return and the current volatility,
that should have mean 0 because sigma squared
t, the t volatility, squared is equal to the square--
or is equal to the expectation of the squared excess residual
return, epsilon t squared.
So if we plug that in, we basically
get an ARMA model for the squared residuals.
And so epsilon t squared is alpha 0 plus alpha 1
plus beta 1 the squared t minus 1 lag plus ut minus beta 1 ut.
And so what this implies is an ARMA 1,1 model with white noise
that has mean 0 and variance 2 sigma to the fourth.
Just plugging things in.
And through our knowledge, understanding,
of univariate time series models, ARMA models,
we can express this ARMA model for the squared residuals
as basically a polynomial lag of the squared residuals
is equal to a polynomial lag of the innovations.
And so we have this expression for what the innovations are.
And it's required that the roots of this a of l operator,
when it thought of it on the complex plane,
have roots outside the unit circle, which corresponds
to alpha 1 plus beta 1 being less than 1 in magnitude.
So in order for these volatility models
not to blow up and be stationary, covariant
stationary, we have these bounds on the parameters.
OK, let's look at the unconditional volatility
or long run variance of the GARCH model.
If you take expectations on both sides of the GARCH model
equation, you basically have the expectation of sigma
squared sub t-- in the long run is sigma star squared-- is
alpha 0 plus alpha 1 plus beta 1 sigma star squared.
So that sigma star squared there is the expectation
of the t minus 1 volatility squared in the limit.
And then you can just solve for this
and see the sigma star squared is equal alpha
0 over 1 minus alpha 1 minus beta 1.
And in terms of the stationarity conditions for the process,
if the long run variance, in order for that to be finite,
you need alpha 1 plus beta 1 to be less than 1 in magnitude.
And if you consider the general GARCH p,1 model,
then the same argument leads to a long run variance being equal
to alpha 0, the sort of intercept term in the GARCH
model, divided by 1 minus the sum of all the parameters.
So these GARCH models lead to constraints
on the parameters that are important to incorporate when
we're doing any estimation of these underlying parameters.
And it does complicate things, actually.
So with maximum likelihood estimation,
the routine for maximum likelihood estimation
is the same for all models.
We basically want to determine the likelihood function
of our data given the unknown parameters.
And the likelihood function is the probability density
function of the data conditional on the parameters.
So our likelihood function as a function
of the unknown parameters c, alpha, and beta
is the value of the probability density, the joint density
of all the data conditional on those parameters.
And that joint density function can
be expressed as the product of successive conditional
expectations of the time series.
And those conditional densities are normal random variables.
So we can just plug in what we know
to be the probability densities of normals
for the t-th innovation epsilon t.
And we just optimize that function.
Now, the challenge with estimating these GARCH models
in part is the constraints on the underlying parameters.
Those need to be enforced.
So we have to have that the alpha i are greater than 0.
Also, the beta j are greater than 0.
And the sum of all of them is between 0 and 1.
Who in this class has had courses in numerical analysis
and done some optimization of functions?
Non-linear functions?
Anybody?
OK.
Well, in addressing this kind of problem, which will come up
with any complex model that you need to estimate,
say via maximum likelihood, the optimization methods
do really well if you're optimizing a convex function,
finding the minimum of a convex function.
And it's always nice to do minimization
over sort of an unconstrained range of underlying parameters.
And one of the tracks in solving these problems
is to transform the parameters to a scale
where they're unlimited in range, basically.
So if you have a positive random variable,
you might use to log of that variable as the thing
to be optimizing over.
If the variable's between 0 and 1,
then you might use that variable divided
by 1 minus that variable and then take the log of that.
And that's unconstrained.
So there are tricks for how you do this optimization, which
come into play.
Anyway, that's the likelihood with the normal distribution.
And we have computer programs that
will solve that directly so we don't
have to worry about this particular case.
Once we fit this model, we want to evaluate
how good it is and the evaluation is based
upon looking at the residuals from the model.
So what we have are these innovations,
epsilon hat t, which should be distributed
with variance or volatility sigma hat t.
Those should be uncorrelated with themselves or at least
to the extent that they can be.
And the squared standardized residuals
should also be uncorrelated.
What we're trying to do with these models
is to capture the dependence, actually,
in this squared residuals, which is measuring
the magnitude of the excess in returns.
So though should be uncorrelated.
There are various test for normality.
I've listed some of those that are the most popular here.
And then there's issues of model selection
for deciding sort of which GARCH model to apply.
I wanted to go through an example of this analysis
with the Euro-dollar exchange rate.
So let me go to this case study note.
So let's see.
There's a package in r called rugarch for univariate GARCH
models, which fits various GARCH models with different--
and fits them by maximum likelihood.
So with this packet-- with this particular library in r,
I fit the GARCH model after actually
fitting the mean process for the exchange rate returns.
Now, when we looked at things last time,
we basically looked at modeling the squared returns.
In fact, there may be an underlying mean process
that needs to be specified as well.
So in this section of the case note,
I initially fit an auto regressive process
using the Akaike information criterion
to choose the order of the auto regressive process
and then fit a GARCH model with normal GARCH terms.
And this is a plot of the normal q,
q plot of the auto regressive residuals.
And what you can see is that the points lie along
a straight line sort of in the middle of the range.
But on the extremes, they depart from that straight line.
This basically is a measure of standardized quantiles.
So in terms of standard units away
from the mean for the residuals, we
tend to get many more high values and many more low values
with the Gaussian distribution.
So that really isn't fitting very well.
If we proceed and fit-- OK, actually that plot
was just the simple ARCH model with no GARCH terms.
And then this is the graph of the q, q plot with the Gaussian
assumption.
So here we can see that the residuals from this model
are suggesting that it may do a pretty good job when things are
only a few standard deviations away from the mean.
Less than two 2, 2.5.
But when we get to more extreme values,
this isn't modeling things well.
So one alternative is to consider a heavier
tailed distribution than the normal, namely
the t distribution.
And consider identifying what t distribution best
fits the data.
So let's just look at what ends up being the maximum likelihood
estimate for the degrees of freedom parameter, which
is 10 degrees of freedom.
This shows the q, q plot when you
have a non-Gaussian distribution that's
t with 10 degrees of freedom.
It basically is explaining these residuals quite well,
so that's accommodating the heavier tailed distribution
of these values.
With this GARCH model, let's see--
if you compare sort of estimate of volatility under the GARCH
and ARCH models-- the GARCH models with the t
distribution-- sorry t distribution versus Gaussian.
Here's just a graph showing time series plots of the estimated
volatility over time, which actually look quite close.
But when you look at the differences,
there really are differences.
And so it turns out that the volatility function
or the volatility estimate GARCH models with Gaussian
versus GARCH with t distributions
are really very, very similar.
And the heavier tailed distribution
of the t distribution means that the distribution
of actual volatility is greater.
But in terms of estimating the volatility,
you have quite similar estimates of the volatility coming out.
And this display-- which you'll be
able to see more clearly in the case notes that I'll post up--
but show that these are really quite similar in magnitude.
And the value at risk concept that was just-- by Ken couple
weeks ago in his lecture from Morgan Stanley--
concerns the issue of estimating what
is the likelihood of returns exceeding some threshold.
And if we use the t distribution for measuring variability
of the excess returns, then the computations in the notes
indicate how you would compute these value at risk limits.
If you compare the t distribution
with a Gaussian distribution at these nominal levels for value
at risk of like 2.5% or 5%, surprisingly you
won't get too much difference.
It's really in looking at sort of the extreme tails
of the distribution that things come into play.
And so I wanted to show you how that plays out
by showing you another graph here.
Those of you who have had a statistics course before
have heard that sort of a t distribution
can be a good approximation to a normal--
or it can be approximated well by a normal
if the degrees of freedom for the t are at some level.
And who wants to suggest a degrees of freedom
that you might have before you're
comfortable approximating a t with a normal?
Danny?
AUDIENCE: 30 or 40.
PROFESSOR: 30 or 40.
Sometimes people say even 25.
Above 25, you can almost expect the t distribution
to be a good approximation to the normal.
Well, this is a graph the PDF for a standard normal
versus a standard t with 30 degrees of freedom.
And you can see that the density functions are very, very close.
The standard-- the CDFs, the Cumulative Distribution
Functions, which is the likelihood
of being less than or equal to the horizontal value
ranges between 0 and 1-- is almost indistinguishable.
But if you look at the tails of the distribution,
here I've computed the log of the CDF function.
You basically have to move much more
than two standard deviations away from the mean
before there's really a difference in the t
distribution with 30 degrees of freedom.
Now I'm going to page up by reducing
the degrees of freedom.
Let's see.
If we could do a page down here.
Page down.
Oh, page up.
OK.
So here is 20 degrees of freedom.
Here's 10 degrees of freedom, in our case, which turns out
to be sort of the best fit of the t distribution.
And what you can see is that, in terms of standard deviation
units, up to about two standard deviations below the mean,
we're basically getting virtually the same probability
mass at the extreme below.
But as we go to four or six standard deviations,
then we get heavier mass with the t distribution.
In discussion of results in finance
when you sort of fit models, people talk about, oh, there
was six standard deviation move or-- which
is just virtually impossible to occur.
Well, with t distributions a six standard deviation move
occurs about 1 in 10,000 times according to this fit.
And so it actually is a common [? idiomatic. ?]
And so it's important to know that these t distributions are
benefiting us by giving us a much better gauge of what
the tail distribution is like.
And we call these distributions leptokurtic,
meaning they're heavier tailed than a normal distribution.
Actually, lepto means slender, I believe,
if you're Greek or have the Greek origin of the word.
And you can see that the blue curve, which
is the t distribution, is sort of a bit more slender
in the center of the distribution, which
allows it to have heavier tails.
All right.
So t distributions are very useful.
Let's go back to this case note here
which discusses-- this case note goes through,
actually, fitting the t distribution--
identifying the degrees of freedom for this t model.
And so with ru GARCH package, we can get the log likelihood
of the data fit under the t distribution assumption.
And here's a graph of the negative log likelihood
versus the degrees of freedom in the t model.
So with maximum likelihood we identify the value,
which minimizes the negative log likelihood.
And that comes out as that 10 value.
All right.
Let's go back to these notes and see
what else we want to talk about.
All right.
OK, with these GARCH models we actually
are able to model volatility clustering.
And volatility clustering is where, over time, you
expect volatility to be high during some periods
and to be low during other periods.
And the GARCH model can accommodate that.
So large volatilities tend to be followed
by large, small volatilities tend
to be followed by small ones.
OK.
The returns have heavier tails than Gaussian distributions.
Actually, even if we have Gaussian errors in the GARCH
model, it's still heavier tailed than a Gaussian.
The homework goes into that a little bit.
And the-- well, actually one of the original papers
by Engle with Bollerslev, who introduced the GARCH model,
discusses these features and how useful
they are for modeling financial time series.
Now, a property of these models that may be obvious, perhaps,
but it is-- OK, these are models that
are appropriate for modeling covariant stationary time
series.
So the volatility measure, which is
a measure of the squared excess return,
is basically a covariant stationary process.
So what does that mean?
That means that's going to have a long term mean.
So with these are GARCH models that are covariant stationary,
there's going to be a long term mean of the GARCH process.
And this discussion here details how this GARCH process
is essentially a mean reversion of the volatility
to that value.
So basically, the sort of excess volatility
of the squared residuals relative to their long term
average is some multiple of the previous period's
excess volatility.
So if we build forecasting models of volatility
with GARCH models, what's going to happen?
Basically, in the long run we predict that any volatility
value is going to revert to this long run average.
And in the short run, it's going to move incrementally
to that value.
So these GARCH models are very good for describing volatility
relative to the long term average.
In terms of their usefulness for prediction,
well, they really predict that volatility
is going to revert back to the mean at some rate.
And the rate at which the volatility reverts back
is given by alpha 1 plus beta 1.
So that number, which is less than 1
for covariant stationarity, is sort of measuring, basically,
how quickly you are reverting back to the mean.
And that sum is actually called a persistence parameter
in GARCH models as well.
So is volatility persisting or not?
Well, the larger alpha 1 plus beta 1 is.
The more persistent volatility is, meaning it's
reverting back to that long run average very, very slowly.
In the implementation of volatility estimates
with the risk metrics methodology,
they actually don't assume that there is a long run volatility.
And so that basically you'll have alpha 1 be equal to 0
and beta 1 equal to, say, 0.95.
So or rather the alpha 0 is 0 and the alpha 1 and beta 1
will actually sum to 1.
And so you actually are tracking a potentially non-stationary
volatility, which allows you to be estimating the volatility
without presuming a long run average is
consistent with the past.
There are many extensions of the GARCH models.
And there's wide literature on that.
For this course, I think it's important to understand
the fundamentals of these models in terms of how they're
specified under Gaussian and t assumptions.
Extending them can be very interesting.
And there are many papers to look at for that.
OK.
let's pause for a minute and get the next topic.
All right.
The next topic is time series, multivariate time series.
In two lectures ago of mine, we talked about univariate time
series and basic methodologies there.
We're now going to be extending that to multivariate
time series.
Turns out there's a multivariate Wold representation theorem,
extension of the univariate one.
There are auto regressive processes
for multivariate cases, which are vector auto regressive
processes.
Least squares estimation comes into play.
And then we'll see where our regression analysis
understanding allows us to specify these vector
auto regressive processes nicely.
There's an optimality properties of ordinary least squares
estimates component wise, which we'll highlight in about a half
an hour.
And go through the maximum likelihood estimation model
selection methods, which are just
very straightforward extensions of the same concepts
for univariate time series and univariate regressions.
So let's talk-- let's introduce the notation
for multivariate time series.
We have a stochastic process, which now is multivariate.
So we have bold x of t is some m dimensional valued
random variable.
And it's stochastic process that varies over time t.
And we can think of this as m different time series
corresponding to the m components of the given
process.
So, say, with exchange rates we could
be modeling m different exchange rate values and want
to model those jointly as a time series.
Or we could have collections of stocks that we're modeling.
So each of the components individually
can be treated as univariate series with univariate methods.
With the multivariate case, we extend
the definition of covariance stationarity
to correspond to finite bounded first and second order moments.
So we need to talk about the first order
moment of the multivariate time series.
Mu now is an m vector, which is the vector of expected values
of the individual components, which we can denote by mu 1
through mu m.
So we basically have m vectors for our mean.
Then for the variance/covariance matrix,
let's define gamma 0 to be the variance/covariance matrix
of the t-th observation of our multivariate process.
So that's equal to the expected value
of xt minus mu xt minus mu prime.
So when we write that down, we have xt minus mu.
This is basically an m by 1 vector
and then xt minus mu prime is a 1 by m vector.
And so the product of that is an m by m quantity.
So the 1, 1 element of that product is the variance of x1t.
And the diagonal entries are the variances
of the components series.
And the off diagonal values are the covariance
between the i-th row series and the j-th column series,
as given by the i-th row of x and the j-th column
of x transpose.
So we're just collecting together all the variances
covariances together.
And the notation is very straightforward and simple
with the matrix notation given here.
Now, the correlation matrix, r0, is obtained by pre and post
multiplying this covariance matrix gamma
0 by a diagonal matrix with the square roots
of the diagonal of this matrix.
Now what's a correlation?
Correlation is the correlation between two random variables
where we've standardize the variables to have
mean 0 and variance 1.
So what we want to do is basically divide through all
of these variables by their standard deviation
and compute the covariance matrix on that new scaling.
That's equivalent to just pre and post multiplying
by that diagonal of the inverse of the standard deviations.
So with matrix algebra, that formula
is-- I think it's very clear.
And this is-- now with-- the previous discussion was just
looking at the sort of contemporaneous covariance
matrix of the time series values at the given time
t with itself.
We want to look at, also, the cross covariance matrices.
So how are the current values of the multivariate time series
xt-- how do they covary with the k-th lag of those values?
So gamma k is looking at how the current period is
vector values as covaried with the k-th lag of those values.
So this covariance matrix has covariance elements
given in this display.
And we can define the cross correlation matrix
by similarly pre and post multiplying
by the inverse of the standard deviations.
The diagonal of gamma 0 is the covariance--
or is the matrix of diagonal entries of variances.
Now, properties of these matrices is-- OK,
gamma 0 is a symmetric matrix that we had before.
But gamma k where k is greater than 1 or less than--
or greater or equal to 1 or less than-- basically different
from 0.
This is not symmetric.
Basically, you may have lags of some variables that
are positively correlated with others and not vice versa.
So the off diagonal entries here aren't necessarily even
of the same sign, let alone equal and symmetric.
So with these covariance matrices,
one can look at how things covary
and whether they are-- whether there is, basically,
a dependence between them.
And you can define-- it's basically
the j star series-- the j star component
of the multivariate time series may
lead the j-th one if the covariance of the k-th lag of j
star is different from 0-- or the covariance of j star k
lags ago is non-zero covaries with the j-th lag.
Sorry.
The current lag.
So xt j star will lead xtj.
Basically, there's information in the lagged values
of j star for the component j.
So if we're trying to build models-- linear regression
models, even, where we're trying to look at how-- trying
to predict values, then if there's a non-zero covariance,
then we can use those variables' information
to actually project what the one variable is given the other.
Now, it can be the case that you have
non-zero covariance in both directions.
And so that suggests that there can
be sort of feedback between these variables.
It's not just that one variable causes another,
but there can actually be feedback.
In economics and finance, there's
a notion of Granger causality.
And basically that-- well, Granger and Engle
got the Nobel Prize number of years ago based on their work.
And that work deals with identifying, in part, judgments
of causality between-- or Granger
causality between variables economic time series.
And so Granger causality basically
is sort of positive or non-zero correlation between variables
where lags of one variable will cause another or cause changes
in another.
All right.
I want to just alert you to the existence of this Wold
decomposition theorem.
This is an advanced theorem, but it's a useful theorem
to know exists.
And this extends-- the univariate Wold decomposition
theorem, which concerns the-- whenever
we have a covariant stationary process,
there exists a representation of that process, which
is the sum of a deterministic process
and s moving average process of a white noise.
So if you're modeling a time series
and you're going to be specifying
a covariant stationary process for that,
there does exist a Wold decomposition representation
of that.
You can basically determine-- identify
the deterministic process that the process might follow.
It might be a linear trend over time or an exponential trend.
And if you remove that sort of deterministic process vt,
then what remains is a process that
can be modeled with a moving average of white noise.
These.
Now here, everything is changed from univariate case
to multivariate case, so we have matrices in place of constants
from before.
So these-- new concepts here are we
have a multivariate white noise process.
That's going to be a process a to t which
is m dimensional which has mean 0.
And the variance matrix of this m vector
is going to be sigma, which is now m by m variance/covariance
matrix of the components.
And that must be a positive semi-definite definite.
And for white noise, we have covariances between, say,
the current t innovation and a lag of its value are 0.
So these are uncorrelated multivariate white noise
processes.
And so they're uncorrelated with each other at various lags.
And the innovation a to t has a covariance of 0
with the deterministic process.
Actually, that's pretty much a given
if we have a deterministic process.
Now, the term psi k-- basically we have this vector xt
is equal to the sum m vectored process vt
plus this weighted average of innovations.
What's required is that the sum of this-- basically each term
psi k and its transpose converges.
Now, if you were to take that xt process
and say let me compute the variance/covariance matrix
of that representation, then you would basically
get terms in the covariance matrix
which includes this sum of terms.
So that sum has to be finite in order for this
to be covariant stationary.
AUDIENCE: [INAUDIBLE].
PROFESSOR: Yes?
AUDIENCE: Could you define what you mean by innovation?
PROFESSOR: Oh, OK.
Well, the innovation is-- let's see.
With-- let me go back up here.
OK.
The innovation process-- innovation process.
OK, if we have, as in this case, we
have sort of our xt stochastic process.
And we have sort of, say, f sub t minus 1
equal to the information on x2 minus 1 x2 minus 2.
Basically consisting of the information set available
before time t.
Then we can model xt to be the expected value of xt given ft
minus 1 plus an innovation.
And so our objective in these models
is to be thinking of how is that process evolving where we can
model the process as well as possible using information up
to time before t.
And then there's some disturbance about that model.
There's something new that's happened
at time t that wasn't available before.
And that's this innovation process.
So this representation with the Wold decomposition
is converting the-- or representing, basically,
the bits of information that are affecting the process that
are occurring at time t and wasn't available prior to that.
All right.
Well, let's move on to vector auto regressive processes.
OK, this representation for a vector auto regressive process
is an extension of the univariate auto regressive
process to p dimensions.
Sorry, to m dimensions.
And so our xt is an m vector.
That's going to be equal to some constant vector c plus a matrix
phi 1 times lag of xt first order, xt minus 1.
Plus another matrix, phi 2 times the second lag of xt,
xt minus 2.
Up to the p-th term, which is a phi p m by m matrix times
x2 minus p plus this innovation term.
So this is essentially-- this is basically how a univariate auto
regressive process extends to an m variate case.
And what this allows one to do is
model how a given component of the multivariate series--
like how one exchange rate varies depending
on how other exchange rates might vary.
Exchange rates tend to co-move together in that example.
So if we look at what this represents
in terms of basically a component series,
we can consider fixing j, a component
of the multivariate process.
It could be the first, the last, or the j-th,
somewhere in the middle.
And that component time series-- like
a fixed exchange rate series or time
series, whatever we're focused on in our modeling--
is a generalization of the auto regressive model
where we have the auto regressive terms
of the j-th series on lags of the j-th series up to order p.
So we have the univariate auto regressive model,
but we also add to that terms corresponding
to the relationship between xj and xj star.
So how does xj, the j-th component,
depend on other variables, other components
of the multivariate series.
And those are given here.
So it's a convenient way to allow for interdependence
among the components and model that.
OK.
This slide deals with representing a p-th order
process as a first order process with vector auto regressions.
Now the concept here is really a very powerful concept that's
applied in time series methods, which
is when you are modeling dependence that goes back, say,
a number of lags like p lags, the structure
can actually be re-expressed as simply a first order
dependence only.
And so it's much easier sort of to deal with just a lag one
dependence then to consider p lag dependence
and the complications involved with that.
So-- and this technique is one where,
in the early days of fitting, like auto regressive
moving average processes and various smoothing methods,
the model-- basically accommodating
p lags complicated the analysis enormously.
But one can actually re-express it just
as a first order lag problem.
So in this case, what one does is one considers
for a vector auto regressive process of order of p,
simply stacking the values of the process.
So let me just highlight what's going on there.
So if we have basically-- OK, so if we have x1,
x2, xn, which are all m by 1 values,
m vectors of the stochastic process.
Then consider defining zt to be equal to xt transpose xt
minus 1 transpose up to x t minus p minus 1 transpose.
Or this is t minus p minus 1.
So there are p terms.
And then if we consider the lagged value of that, that's
x2 minus 1, x2 minus 2, x2 minus p transpose.
So what we've done is we're considering zt.
This is going to be m times p.
It's actually 1 by m times p in this notation.
Well, actually I guess I should put transpose here.
So m minus p by 1.
OK, in the lecture notes it actually
is primed there to indicate the transpose.
Well, if you define zt and zt minus 1 this way,
then zt is equal to d plus a of zt minus 1
plus f, where this is d.
Basically the constant term has the c entering and then 0's
everywhere else.
And the a matrix is phi 1 phi 2 up to phi p.
And so basically the zt vector transforms the zt--
or is the transpose-- this linear transformation
of the zt minus 1.
And we have sort of a very simple form
for the constant term and a very simple form for the f vector.
And this is-- renders the model into a sort of a first order
time series model with a larger multivariate series,
basically mp by 1.
Now, with this representation we basically have-- we
can demonstrate that the process is going to be stationary
if all eigenvalues of the companion matrix a
have modulus less than 1.
And let's see-- if we go back to the expression.
OK, if the eigenvalues of this matrix A are less than 1,
then we won't get sort of an explosive behavior
of the process when this basically increments over time
with every previous value getting multiplied by the A
matrix and scaling the process over time by the A-th power.
So that is required.
All eigenvalues of A have to be less than 1.
And equivalently, all roots of this equation
need to be outside the unit circle.
You remember there was a constraint of-- or a condition
for univariate auto regressive models
to be stationary, that the roots of the characteristic equation
are all outside the unit circle.
And the class notes go through and went
through the derivation of that.
This is the extension of that to the multivariate case.
And so basically one needs to solve
for roots of a polynomial in Z and determine
whether those are outside the unit circle.
Who can tell me what the order of the polynomial
is here for this sort of determinant equation?
AUDIENCE: [INAUDIBLE] mp.
PROFESSOR: mp.
Yes.
It's basically of power mp.
So in a determinant you basically
are taking products of the m components
in the matrix, various linear combinations of those.
So that's going to be an mp dimensional polynomial.
All right.
Well, the mean of the stationary VAR process
can be computed rather easily by taking expectations
of this on both sides.
So if we take the expectation of xt
and take expectations across both sides,
we get that mu is the c vector plus the product of the phi
case times mu plus 0.
So mu, the unconditional mean of the process,
actually has this formula just solving
for mu in the top-- in the second line to the third line.
So here we can see that basically this expression
1 minus phi 1 through phi p, that inverse has to exist.
And actually, if we then plug in the value
of c in terms of the unconditional mean,
we get this expression for the original process.
So the unconditional mean c, if we demeaned the process,
there's busy know mean term.
There's 0.
And so basically the mean adjusted process x
follows this multivariate vector auto regression
with no mean, which is actually used when this is specified.
Now, this vector auto regression model
can be expressed as a system of regression equations.
And so what we have with the multivariate series, if we have
multivariate data, we'll have n sample observations,
xt which is basically the m vector
of the multivariate process observed for n time points.
And for the computations here, we're
going to assume that we have p sort of-- we
have pre-sample observations available to us.
So we're essentially going to be considering models
where we condition on the first p time
points in order to facilitate the estimation methodology.
Then we can set up m regression models corresponding
to each component of the m variate series.
And so what we have is our original--
we have our collection of data values, which is x1 transpose,
x2 transpose, down to xn transpose,
which is an n by m matrix.
OK, this is our multivariate time series
where we were just-- the first row corresponds
to the first time values, nth row to the nth time values.
And we can set up m regression models
where we're going to consider modeling
the j-th column of this matrix.
So we're just picking out the univariate time series
corresponding to the j-th component.
That's yj.
And we're going to model that has z beta j plus epsilon j
where z is given by the vector of lagged values
of the multivariate process where there's,
for the t-th-- t minus first value
we have that current value-- or the t
minus first, t minus second, up to t minus p.
So we have basically p m vectors here.
And so this j-th time series has elements
that follow a linear regression model
on the lags of the entire multivariate series up to p
lags with their progression parameter given by beta j.
And basically the beta j regression parameters
corresponds to the various elements of the phi matrices.
So now there's a one-to-one one correspondence between those.
All right.
So I'm using now a notation where superscript j corresponds
to the j-th component of the series,
of the multivariate stochastic process.
So we have an mp plus 1 vector progression parameters
for each series j, and we have an epsilon j
for an n-vector innovation errors for each series.
And so basically if this, the j-th column, is yj,
we're modeling that to be equal to the simple matrix
Z times beta j plus epsilon j, where this is n by 1.
This is n by np plus 1.
And this beta j is the mp plus 1 progression parameter.
OK.
One might think, OK, one can consider
each of these regressions for each of the component series,
you could consider them separately.
But to consider them all together,
we can define the multivariate regression model,
which has the following form.
We basically have the n vectors for the first component,
and then the second component up to nth component.
So an n by p matrix of dependent variables,
where each column corresponds to a different component series,
follows a linear regression model
with the same Z matrix with different regression
coefficient parameters, beta 1 through beta m corresponding
to the different components of the multivariate series.
And we have epsilon 1, epsilon 2, up to epsilon m.
So we're thinking of taking-- so basically the y1 y2 up to ym
is essentially this original matrix of our multivariate time
series because it's the first component in the first column
and the nth component in the nth column.
And the-- this regression parameter
or this explanatory variables matrix X, Z in this case
corresponds to lags of the whole process up to p lags.
So we're having lags of all the m variate process up to p lags.
So that's mp and then plus 1 for our constant.
So this is the set up for a multivariate regression model.
In terms of how one specifies this,
well, actually, in economic theory
this is also related to seemingly unrelated
regressions, which you'll find in econometrics.
If we want to specify this multivariate model, well,
what we could do is we could actually
specify each of the component models
separately because we basically have sort of-- can think
of the univariate regression model for each component
series.
And this slide indicates basically what
the formulas are for that.
So if we don't know anything about multivariate regression
we can say, well, let's start by just
doing the univariate regression of each component series
on the lags.
And so we get our beta hat j's least squares
estimates given by the usual formula where
the independent variable is matrix Z goes Z transpose Z
inverse Z transpose Y of the residual.
So these are familiar formulas.
And if we did this for each of the component series j,
then we would actually get sample estimates
of the innovation process, eta 1.
Basically the whole eta series.
And we could actually define from these estimates
of the innovations our covariance matrix
for the innovations as the sample covariance
matrix of these etas.
So all of these formulas are-- you're basically
applying very straightforward estimation
methods for the parameters of a linear regression
and then estimating variances/covariances
of these innovation terms.
So from this, we actually have estimates
of this process in terms of the sigma and the beta hats.
But it's made assuming that we can
treat each of these component regressions separately.
A rather remarkable result is that
these component-wise regressions are actually
the optimal estimates for the multivariate regression
as well.
And as mathematicians, this kind of result
is, I think, rather neat and elegant.
And maybe some of you will think this is very obvious,
but it actually-- it isn't quite obvious.
That said, this component-wise estimation
should be optimal as well.
And the next section of the lecture notes
goes through this argument.
And I'm going to, in the interest of time,
go through this-- just sort of highlight what the results are.
The details are in these notes that you can go through.
And I will be happy to go into more detail about them
during office hours.
But if we're fitting a vector auto regression model where
there are no constraints on the coefficient matrices,
phi 1 through phi p, then these component-wise estimates,
accounting for arbitrary covariance matrix sigma
for the innovations, those basically
are equal to the generalized least squares estimates
of these underlying parameters.
You'll recall we talked about the Gauss Markov theorem
where we were able to extend the assumption of sort
equal variances across observations
to unequal variances and covariances.
Well, it turns out to these component-wise OLS
estimates are, in fact, the generalized least squared
estimates.
And under the assumption of Gaussian distributions
for the innovations, they, in fact,
are maximum likelihood estimates.
And this theory applies Kronecker products.
We're not going to have any homework with Kronecker
products.
These notes really are for those who
have some more extensive background in linear algebra.
But it's a very nice use of these Kronecker product
operators.
Basically, this notation-- or no, x circle--
I'll call it Kronecker-- is one where
you take a matrix A and a matrix B
and you consider the matrix which
takes each element of A times the whole matrix B.
So we start with an m by n matrix A
and end up with an mp by qn matrix
by taking each element of A times the whole matrix B.
So it's, they say, has this block structure.
So this is very simple definition.
If you look at properties of transposition of matrices,
you can prove these results.
These are properties of the Kronecker product.
And there's a vec operator which takes a matrix
and simply stacks the columns together.
And in the talk last Tuesday of Ivan's, talking about modeling
the volatility surface, he basically, he
was modeling a two dimensional surface-- or a surface
in three dimensions, but there was
two dimensions explaining it.
You basically can stack columns of the matrix
and be modeling a vector instead of a matrix of values.
So the vectorizing operator allows us to manipulate terms
into a more convenient form.
And this multivariate regression model
is one where it's set up as sort of a n by m matrix Y,
having that structure.
It can be expressed in terms of the linear regression form
as y star equaling the vector, the vec of y.
So we basically have y1 y2 down to ym all lined up.
So this is pm by 1.
That's going to be equal to some matrix plus the epsilon 1
epsilon 2 down to epsilon n.
And then there's going to be a matrix
and a regression coefficient matrix beta 1 beta
2 down to beta p.
So we consider vectorizing the beta matrix,
vectorizing epsilon, and vectorizing y.
And then in order to define this sort
of simple linear regression model, univariate regression
model, well, we need to have a Z in the first column
here corresponding to beta 1 for Y1 and 0's everywhere else.
In the second block we want to have
a Z in the second off diagonal with 0's everywhere else and so
forth.
So this is just re-expressing everything in this notation.
But the notation is very nice because, at the end of the day
we basically have a regression model like we
had when we were doing our regression analysis.
So all the theory we have for specifying these models
plays through with univariate regression.
And one can go through this technical argument
to show that the generalized least squares estimate is,
in fact, the equivalent to the component-wise values.
And that's very, very good.
Maximum likelihood estimation with these models.
Well, we actually use this vectorized notation
to define the likelihood function.
And if these assumptions are made
about the linear regression model,
we basically have an n times m vector
of dependent variable values, whereas your multivariate
normal with mean given by x star beta star
and then a covariance matrix epsilon.
The covariance matrix of epsilon star is sigma star.
Well, sigma star is In Kronecker product sigma.
So if you go through the math of this,
everything matches up in terms of what the assumptions are.
And the conditional probability density function of this data
is the usual functions of log normal or of a normal sample.
So we have unknown parameters beta star sigma,
which are equal to the joint density
of this normal linear regression model.
So this corresponds to what we had
before in our regression analysis.
We just had this more complicated definition
of the independent variables matrix X star.
And a more complicated definition
of our variance/covariance matrix sigma star.
But the log likelihood function ends up
being equal to a term proportional
to the log of the determinant of our sigma matrix
and minus one half q of beta sigma, where q of beta sigma
is the least squares criterion for each of the component
models summed up.
So the component-wise maximum likelihood estimation
is-- for the underlying parameters,
is the same as the large one.
And in terms of estimating the covariance matrix,
there's a notion called the concentrated log likelihood,
which comes into play in models with many parameters.
In this model, we have unknown parameters--
our regression parameters beta and our covariance matrix
for the innovations sigma.
It turns out that our estimate of the regression parameter
beta is independent, doesn't depend-- not statistically
independent-- but does not depend
on the value of the covariance matrix sigma.
So whatever sigma is, we have the same maximum likelihood
estimate for the betas.
So we can consider the log likelihood
setting the beta parameter equal to its maximum likelihood
estimate.
And then we have a function that just
depends on the data and the unknown parameter sigma.
So that's a concentrated likelihood function
that needs to be maximized.
And the maximization of the log of a determinant of a matrix
minus n over 2, the trace of that matrix times an estimate
of it, that has been solved.
It's a bit involved.
But if you're interested in the mathematics for how that's
actually solved and how you take derivatives of determinants
and so forth, there's a paper by Anderson and Olkin
that goes through all the details of that
that you can Google on the web.
Finally, let's see.
There's-- well, not finally.
There's model selection criteria that can be applied.
These have been applied before for regression models
for univariate time series model, the Akaike Information
Criterion, the Bayes Information Criterion, Hannan-Quinn
Criterion.
These definitions are all consistent
with the other definitions.
They basically take the likelihood function
and you try to maximize that plus a penalty
for the number of unknown parameters.
And that's given here.
OK, then the last section goes through
an asymptotic distribution of least squares estimates.
And I'll let you read that on your own.
Let's see.
For this lecture I put together an example of fitting vector
auto regressions with some macroeconomic variables.
And I just wanted to point that out to you.
So let me go to this document here.
What have we got here?
All right.
Well, OK.
Modeling macroeconomic time series is an important topic.
It's what sort of central bankers do.
They want to understand what factors are affecting
the economy in terms of growth, inflation, unemployment.
And what's the impact of interest rate policies.
There are some really important papers
by Robert Lederman and Christopher Sims dealing
with fitting vector auto regression
models to a macroeconomic time series.
And actually, the framework within which
they specified these models was a Bayesian framework,
which is an extension of the maximum likelihood method where
you'll incorporate reasonable sort
prior assumptions about what the parameters ought to be.
But in this note, I sort of basically
go through collecting various macroeconomic variables
directly off the web using the package r.
All this stuff is-- these are data
that you can get your hands on.
Here's the unemployment rate from January 1946
up through this past month.
Anyone can see how that's varied between much less than 4%
to over 10%, as it was recently.
And there's also the Fed funds rate,
which is one of the key variables
that the Federal Reserve Open Market Committee controls,
or I should say controlled in the past,
to try and affect the economy.
Now that value of that rate is set almost at zero
and other means are applied to have an impact
on economic growth the economic situation
of the market-- of the economy, rather.
Let's see.
There's also-- anyway, a bunch of other variables.
CPI, which is a measure of inflation.
What this note goes through is the specification
of vector auto regression models for these series.
And I use just a small set of cases.
I look at unemployment rate, federal funds,
and the CPI, which is a measure of inflation.
And there's-- if one goes through,
there are multivariate versions of the autocorrelation
function, as given on the top right panel here,
between these variables.
And one can also do the partial autocorrelation function.
You'll recall that autocorrelation functions
and partial autocorrelation functions
are related to what kind of-- or help us understand what kind
of order ARMA processes might be appropriate
for univariate series.
For multivariate series, then there are basically
cross lags between variables that are important,
and these can call be captured with vector auto regression
models.
So this goes through and shows how
these things are correlated with themselves.
And let's see.
At the end of this note, there are some impulse
response functions graphed, which
are looking at what is the impact of an innovation in one
of the components of the multivariate time series.
So like if Fed funds were to be increased by a certain value,
what would the likely impact be on the unemployment rate?
Or on GNP?
Basically, the production level of the economy.
And this looks at-- let's see.
Well, actually here we're looking
at the impulse function.
You can look at the impulse function of innovations
on any of the component variables on all the others.
And in this case, on the left panel
here is-- it shows what happens when unemployment
has a spike up, or unit spike.
A unit impulse up.
Well, this second panel shows what's
likely to happen to the Fed funds rate.
It turns out that's likely to go down.
And that sort of is indicating-- it's sort
of reflecting what, historically,
was the policy of the Fed to basically reduce interest
rates if unemployment was rising.
And then-- so anyway, these impulse response functions
correspond to essentially those innovation
terms on the Wold decomposition.
And why are these important?
Well, this indicates a connection, basically,
between that sort of moving average representation
and these time series models.
And the way these graphs are generated
is by essentially finding the Wold decomposition
and then incorporating that into these values.
So-- OK, we'll finish there for today.