Description
ECE478 Financial Signal Processing Problem Set I
1. Be a Financial Engineer!
Here we will study mixes of European calls and puts to craft di§erent types of payouts.
Here assume maturity date T, with underlying security price ST at T. The payout for
a European call and put, respectively, at time T with strike price K is, respectively:
c (T; T) = (ST K)
+
p (T; T) = (K ST )
+
DeÖne a digital call or put option to have respective payout as follows:
dC (T; T) = 1 if ST > K, 0 otherwise
dp (T; T) = 1 if ST < K 0 otherwise
A long position in one unit of a derivative with payout DT means the value of the asset
in the portfolio is DT . A short position in one unit means the value is DT .
Combining various amounts of long and short positions in call and put options can yield
a general piecewise linear but overall continuous value V (ST ) at the time of maturity.
Including digital options allows discontinuities. Here we will see several important
combinations.
Use a computer to generate plots by assigning reasonable values to the parameters K,
etc. You just need to produce one example of each. In all cases, horizontal axis is
ST > 0 and the vertical axis is V (ST ).
Note: As you graph these, keep in mind you are looking at the ìpayoutî V (ST ). If
you purchase one of these derivatives at a price say Vt at t, then your net ìproÖtîor
ìloss” at time T would be V (ST ) Vt=Z (t; T) (i.e., subtract o§ the money you would
have had in your pocket if you never purchased the derivative). This would be a more
useful graph, perhaps, e.g., to determine the range of ST for which you end up with a
net positive, but that would require determining the price at t, Vt
, which we are not
dealing with here. One thing to look for in any of these instruments: is the potential
loss unbounded? In other words, is min V (ST ) Önite?
Also, in what follows below, except for one case Iím not asking you to derive an analytic
formula for V (ST ). If you write core code to compute the payout for the call, put,
digital call, digital put, then just write code to combine these as suggested to obtain
the graphs. Your code should be general purpose, e.g., with parameters such as K
adjustable. When you generate graphs, plug in reasonable values, e.g., K = 1 or
whatever. Scale your horizontal and vertical axes accordingly.
(a) Graph V (ST ) for each of the following: long call, long put, short call, short put
(b) A straddle is a combination of one unit of a call and one unit of a put for the
same security at the same strike price K. Verify V (ST ) = jST Kj, and graph
this for one choice of K.
(c) A call-put spread is to long a call at K1 and short a call at K2 with K1 < K2.
Graph V (ST ) for one choice of the pair K1; K2.
(d) A butteráy is a combination of the following. Let K1 < K2, and 0 < < 1. Let
K = K1 + (1 ) K2. Then:
long calls at strike price K1
long (1 ) calls at strike price K2.
short 1 call at strike price K
.
Fix one choice for K1 < K2. Graph a butteráy V (ST ) for each of the following
cases: = 1=3; 1=2, 2=3.
(e) Call ladder: with K1 < K2 < K3, long one K1 call, short one K2 call, short one
K3 call.
(f) A digital call spread: with K1 < K2, a long a digital call with strike K1 and short
a digital call with strike K2.
ECE478 Problem Set II: Portfolio Analysis
Data
Our primary source of data in this experiment will be S&P 500, USD LIBOR rate (which
we will take as the risk-free rate), and stock data you can get from online sources such as
Yahoo Önance.
Select 5 stocks that you expect to have some correlation, e.g., from a common sector. [Make
sure these companies have been publicly traded since 2000]
Take data over the period 2000 through 2018. We will keep data over each year separate.
For purposes of computing the risk-free return, take the average LIBOR over the calendar
year. Note that the deÖnition of the daily return is the annual return (i.e., for a one year
time horizon) divided by 360; this is the daycount conversion for interest rates, not stocks!
If you look at monthly returns, then divide by 12.
Pull daily date for S&P500 and your 5 stocks. Use the adjusted closing price. This compensates for e§ects such as dividend payments.
Compute the daily log returns.
Compute the monthly log returns by comparing the Örst (trading) day of a month to the
Örst (trading) day of the prior month.
For the Örst day (month) of a calendar year, use the Önal day (month) of the previous
calendar year. (E.g., you will need some data from 1999 to do this for the year 2000).
In what follows, you can use ìcannedîroutines to compute mean vectors, covariance matrices,
and linear algebra operations such as eigenanalysis. However, the rest you should build up
yourself. Donít use higher level economic/Önancial functions.
There are potentially many graphs generated. You have to decide how to organize things
intelligently. I also want to see comments, especially when data doesnít seem to match
theory.
Analysis
Write code to take the daily data over a period of one calendar year. Consider a portfolio
comprised of your 5 stocks plus the risk-free return derived from USD LIBOR.
1
When you do this, there may be years when the standard theory fails because the USD
LIBOR is too high. Even if it is close to the assumed threshold, though technically below
the level, results can look weird. Have your code identify those conditions, as an exception.
If the exception occurs, then just pick an artiÖcial but convenient value for the risk-free
return, just to push through the theory. But you should clearly indicate when that happens.
You should perform the following operations for several years (no, not for EVERY year).
1. Compute the mean vector m and covariance matrix C for the stocks.
2. Compute the weight vectors of the MVP and Market portfolios ( ~wMV P , ~wM, respectively), and the respective points (MV P ; MV P ) and (M; M).
3. Generate the points on the e¢ cient frontier. Graph it, and also place markers at the
MVP and market portfolio points. Also superimpose a graph of the CAPM.
4. Identify the portion of the e¢ cient frontier that involves no short selling at all, and
highlight it with a di§erent color.
5. Look at the MVP and Market Portfolio weight vectors. Some of them may be negative.
If this is the case, zero those terms out, and adjust the other weights (rescale them)
to sum to 1. We can take this as an approximation to the optimal solution subject to
the constraint of no short-selling. Place markers for these ìmodiÖedîMVP and MP
points.
6. Pick two stocks out of 5 and compute the MP for just this pair. Compute the correlation
coe¢ cient between this portfolio and the MP (for the full set of 5 stocks), the , the
systematic risk and diversiÖable risk.
7. Now study the covariance matrix C:
(a) Compute the eigenvalues, and show that are all positive.
(b) Find the condition number, max=min. This is one way to quantify how numerically ìunstableîmanipulating C is, e.g., how close to nonsingular it is (the larger
this is, the worse for us).
(c) Graph the eigenvalues on a log scale in descending order.
8. Now we look at the S&P 500. Compute its (; ) and place a marker for it on the
above graph, to compare it to the other portfolios you have been looking at.
9. Compare the S&P500 with the MP for the full 5 stocks, and the MP for just the 2
stock subset you took, and also for the ímodiÖed” MP in which you removed any short
selling. . In particular, do one of these dominate some of the others? Compute the
correlation coe¢ cient between the S&P500 and these others.
2
ECE478 Problem Set III: Binomial Asset Pricing Model
Write code to simulate the BAPM and perform analysis as suggested. Notation and so forth
follows as given in Shreve, Stochastic Calculus for Finance vol. 1. Your simulation should
take r; d; u as given, and constant at all time steps. Building in error detection to ensure
these satisfy the no-arbitrage constraint is optional. Might as well take S0 = 1. We will
assume the coin toss distribution is constant over time, but the particular (p; q) values can
be variable. Two special cases are the ìactualî probabilities (p0; q0), and the risk-neutral
measure (~p; q~). Let Fn be the -algebra generated by the Örst n tosses, which covers the
time span f0 k ng. Let ~!n = (!1 !n) 2 fH; Tg
n
denote a particular ëpathíthrough
the Örst n tosses. Assume every stochastic process here, say Xn is adapted to this Öltration.
Your model should cover the time steps f0 n Ng. When you perform a Monte Carlo
simulation, M is the number of samples used (i.e., for purposes of Monte Carlo, you repeat
the experiment M times and average over those results).
There are two types of derivatives: one type has a payout VN = VN (SN ) that depends only
on the Önal value of the security; this is called ëpath independentí. Another type can depend
on the values of the security over the whole time span, i.e., VN = VN (S0; S1; ; SN ). Here
we will only build out simulation for the case of path-independent payouts, so feel free to
write your code assuming that condition holds. [A more complete BAPM simulation would
have to be able to handle the more general case.]
Finally, in what follows, for any stochastic process Xn, the discounted process is X~
n =
1
(1+r)
n Xn.
1. Distributions: Before you start coding: Assuming the coin-toss distribution (p; q),
Önd the distribution of Sn. That is, list all the possible values (there are n+1 of them),
and their probabilities. Also, if Rn = log (Sn=S0), show that Rn = cnYn + dn where
Yn binomial (n; p) (Önd the constants cn; dn). Note: In the continuous-time case,
something similar will occur: we will get Rt = log (St=S0) Gaussian (what else?),
and hence the distribution of St will be what is called lognormal.
2. Exact simulation: Assume we have code to compute the payout function V (SN ) for
a path-independent derivative.
(a) Write code to compute V0 = Ep~
V~N
using the known distribution for SN .
(b) Write code to generate one step of the replicating portfolio. If at time n there
are n shares of the stock, and wealth Xn, then the amount Mn = Xn nSn is
held in the money market. From n ! n + 1, the stock price changes from Sn to
Sn+1, and the amount in the money market grows to (1 + r) Mn. Thus, at time
n + 1, the wealth equation states:
Xn+1 = nSn+1 + (1 + r) (Xn nSn)
For this to be a replicating portfolio, Xn = Vn, the price of the derivative at time n,
where at n = N this should matchVN (SN ) exactly. Give ~!n, 0 n N 1, there
are two possibities for ~!n+1, namely (~!nH) and (~!nT). Assuming Vn+1 (~!nH),
Vn+1 (~!nT) are both known, and of course Sn+1 (~!nH); Sn+1 (~!nT); Sn (~!n) are
all known, with given r, the wealth equation yields two linear equations in the
unknowns n (~!n); Xn = Vn (~!n). By running this recursively backwards, you
should be able to derive fn (~!n)gall ~!n, for 0nN1
, and of course X0 = V0. Write
code to do this! Also as you do this, keep track of whether any n or Xn nSn
values are negative (the Örst is short selling the stock, the second is ìshort selling
the money marketî, i.e., borrowing money to buy stock).
(c) Now let us take a speciÖc example to test your code. Take r = 0:05, u = 1:1,
d = 1:01, N = 5, and the derivative a European call option with strike price
K = Ep~ (SN ) (you know the distribution of SN so should be able to Önd a simple
formula for this!) There are only 2
N = 32 paths so it should not be unreasonable
to compute and store the complete set of replicating portfolio values fn (~!n)g.
In any case, compute V0 using each of your two methods, and check that they
match! Question: Do you have to do any short selling the stock or borrowing
from the money market along the way?
3. Monte Carlo simulation: Next is to deal with the case N = 100. There are 2
100
paths, so rather than computing over all paths, we will use a Monte Carlo approach.
So let me Örst describe the code you should create. For Öxed m, and an adapted
stochastic process Xn, we want to compute Yn = Ep~ (Xn+mjFn). What this means
is Yn = Yn (~!n), a function of the path associated with the Örst n tosses. So assume
that is given. You will generate M random paths (!n+1 !n+m) according to the
distribution, compute Xn+m (~!n!n+1 !n+m) for each, and average over them to get
Yn. The details of computing Xn+m should be contained in a function that is called by
the Monte Carlo ìwrapperî. It is ok if you have to modify your wrapper code a little
bit in terms of the variables you pass to or receive back from this core function.
(a) As a Örst step, let us check the underlying behavior of your Monte Carlo code. Go
back to the previous case, with N = 5 and r; u; d as given. Taking M = 1; 5; 10; 32,
estimate S0 from SN , and V0 from VN . Note that you are not going through every
path systematically! You are generating paths independently, and so it is possible
(especially in the low order case) the same path will occur multiple times, so even
with M = 32 there is no guarantee you will get all paths! Hence your S0; V0
estimate will not be exact. The idea is to see how close your estimates are for
various M.
(b) Now take the case N = 100. In this case, we should change some of the parameters. Take r = 103
(so over 100 steps there is a total return of about 10%),
u = 1 + 5 103
, d = 1 + 104
. Take di§erent values of M, starting small so your
computer doesnít take too long to run, and increasing it somewhat. We want to
hopefully see an e§ect of increasing M, without forcing you to run simulations
overnight! First, as a check, we know we should get S0 = Ep~
S~N
and fortunately we know S0 = 1. See how this works for various M. Next, again let V (S
be a European call with K = Ep~ (SN ) [Note: Careful, this is not the discounted
value, since the strike price applies at tine N; also we are computing the expectation assuming p~; there is nothing formal here, this is just a SUGGESTION I am
making so you have a reasonable value K that lies somewhere between min SN
and max SN ] Compute V0 = Ep~
V~N
for the same M as before. Also compute
V10 assuming the Örst ten tosses are all heads, and again assuming all are tails.
Last thing to do: instead of using p~, take two cases for the ìactualîprobabilities:
p0 = 0:9~p, and p0 = 1:1~p. For each case, compute Ep0
S~N
, and compare with
S0, and Ep0
V~N
. In the latter case, if the answer with p0 is below V0 then
you should not have purchased the option (putting the cash into a money market
would have been better), if above V0 then purchasing the option gave you a higher
return than stashing you money in the money market 😀 So last question: even
if I told you that your expected return was better than putting it into a money
market, you may still have decided not to purchase this option. Why not? Hint:
Even if you know p0 (which in this case would lead to this result), we still donít
have an arbitrage here. What does that mean?
3
ECE478 Problem Set IV: Stochastic Calculus
Notation here follows Shreve, Stochastic Calculus for Finance, vols. I & II, Springer, 2004.
Theoretical Problems
1. Let Ft be the Öltration generated by a Wiener process W (t). Let R (t) be the interest
rate process used to deÖne the discount process D (t). Assume there exists a unique
risk-neutral measure, leading to the Wiener process W~ (t) with respect to P~. If V (T)
is a random variable that is FT -measurable, and V (t) is deÖned via:
V (t) = 1
D (t)
E~ (D (T) V (T)jFt)
then D (t) V (t) is a martingale.
(a) Suppose V (T) > 0 a.s. Show that V (t) > 0 a.s. (from its deÖnition above).
(b) Show that there exists an adapted process ( ~ t) such that:
dV (t) = R (t) V (t) dt +
( ~ t)
D (t)
dW~ (t)
Hint: Start with a formula for d (D (t) V (t)) as per the martingale representation
theorem, then as you expand this out recognize that dV dt = 0.
(c) Show that there exists an adapted process (t) such that we can write:
dV (t) = R (t) V (t) dt + (t) V (t) dW~ (t)
By the way, (t) can be random and in particular it is Öne if the formula for
(t) you derived involves V (t). The point is there is SOME process you can
put there that works! How did we use strict positivity? (Think of D (t), V (t)
as continuous processes; D (t) is intrinsically positive, but what happens if V (t)
can take on negative as well as positive values?) This shows that V (t) is a
generalized geometric Brownian motion process. The point of this problem is
that every strictly positive asset is a generalized geometric Brownian motion.
2. Let X (t); Y (t) be ItÙ processes given by:
dX (t) = a (t) dt + b (t) dW (t)
dY (t) = c (t) dt + d (t) dW (t)
where a; b; c; d are adapted processes. Assume Y (t) > 0 a.s., and V (t) = X (t) =Y (t).
Obtain an SDE satisÖed by V (t), simpliÖed so it has the above form (i.e., in the form
of an ItÙ process). Note that X (t); Y (t), but not dX (t) or dY (t), can appear in your
Önal expression for dV (t).
Simulating Stochastic Di§erential Equations
We are going to be viewing continuous-time deterministic functions and stochastic processes
in discretized time. For simplicity, we will take equally spaced time intervals, say t = n,
0 n N, with Önal time T = N. (Careful, including t = 0, this is N + 1 points). We
will use the notation:
x [n] = x (n) = xn
Consider a stochastic di§erential equation of the form:
dX (t) = (t; X (t)) dt + (t; X (t)) dW (t)
where (); () are deterministic functions, and as usual W (t) is a Wiener process. The
parameter controls the variance of the increment dW; speciÖcally, in the discretization, let
us use the notation:
dW [n] = W [n + 1] W [n]
where fdW [n]g are iid N (0; ). The SDE actually means:
X (t) = X (0) + Z t
0
(u; X (u)) du +
Z t
0
(u; X (u)) dW (u)
In discretized form, with t = n this becomes:
X [n] = X [0] + Xn1
m=0
(m; X [m]) +
Xn1
m=0
(m; X [m]) dW [m]
Letís default with = 0:01 and N = 250 (which is approximately the number of days per
year for a typical Önancial instrument).
1. Start with basic geometric Brownian motion:
dSt = Stdt + StdWt
with S (0) = 1, = 0:1, = 0:2. Also assume a constant underlying interest rate
r = 0:05. Under the risk-neutral measure, dSt satisÖes a modiÖed SDE involving dt
and dW~
t
. See Problem 3 below: If you ever detect St 0 trap that condition,
note that it occurred, and discard that path (donít use it).
(a) Write the modiÖed SDE involving dW~
t
. SpeciÖcally, calculate the coe¢ cients that
appear in this SDE from the speciÖc values of ; ; r provided. Unless speciÖed
otherwise, however, grow your paths using the ORIGINAL formulation (i.e., the
ACTUAL probabilities).
(b) Generate 1000 paths of S [n].
(c) Use a Monte Carlo approach to estimate E (S [N=2]) and E (S [N]) directly.
(d) We now want to connect this to the Black-Scholes model. Let V (T) = (S (T) K)
+
,
the payout of a European call option with strike price K. In our discrete notation,
V [N] = (S [N] K)
+
. Use your estimate for E (S [N]) as your value for K. Use
the Black-Scholes model to obtain a formula for V [N=2] in terms of S [N=2], and
graph it (for the speciÖed parameters ; ; r
(e) Now take the Örst 10 paths you generated, S
(i)
[n], 1 i 10. For each S
(i)
[N=2],
you can compute V
(i)
[N=2] from the Black-Scholes formula. On the other hand,
you can use a Monte Carlo approach to compute V
(i)
[N=2] using the martingale
property of the discounted stock price. So, for each i, grow 1000 paths from N=2
to N and average to estimate V
(i)
[N=2]. Compare these estimated values with the
exact values in these cases. Report the results how you see Öt: a table; you could
superimpose two scatter plots (S
(i)
[N=2] versus actual V
(i)
[N=2], and S
(i)
[N=2]
versus actual V
(i)
[N=2]). Comments: Careful. First, you need to grow NEW
paths, emanating from time t = T=2, towards t = T, taking S
(i)
[N=2] as an initial
condition. Also, you are computing V
(i)
[N=2] as an expectation by virtue of the
martingale property, BUT the question is WHICH SDE FOR St DO YOU
USE?
2. Cox-Ingersoll-Ross Interest Rate Model
dR (t) = (a R (t)) dt +
p
R (t)dW (t)
with ; ; positive, and let R (0) = r > 0. This model for interest rate R (t) guarantees R (t) > 0. This model is typically used for short term interest rates, which tend to
exhibit volatility. Although no closed form solution exists, it is possible to determine
certain properties for it. In particular the mean and variance of R (t) are given by:
E (R (t)) = e
tr +
1 e
t
!
var (R (t)) =
2
r
e
t e
2 t
+
2
2
2
1 2e
t e
2 t
!
2
2
2
Discretizing this SDE for simulation purposes is very tricky because the property that
R (t) > 0 does not carry over if you use the simple discretization as in the previous problem! This actually brings up a larger questionñ does the discretized system
accurately reáect properties of the original SDE?
(a) Select: = 1, = 0:10 , r = 0:05, = 0:5. Generate 1000 paths for R (t) over
a span 0 t 10, using = 0:01. In your code, trap the condition that R 0
(if this occurs, your code should display an exception, should halt that one path,
but should continue computing the other paths). So the goal is to generate 1000
valid paths (with R (t) never reaching 0). How many ìbadî paths do you get
in order to reach 1000?
(b) Graph the Örst 10 paths for R (t), just to see what it looks like.
(c) Use a Monte Carlo approach to estimate the mean and variance of R (t) at t = 1
and t = 10, and compare with the exact formulas given above. Here you should
use the valid paths only.
3. Go back to Problem 1 again. We know the actual St (in continuous-time) is geometric
Brownian motion so we necessarily have St > 0 a.s., and yet is that property guaranteed in your simulation? On the other hand, with all your simulations, did you ever
encounter simulated St 0? I expect not. Why not? [Even if you got a negative
answer, why did I guess that you wouldnít? The hint is that it is a guess
ECE478 Problem Set V: Financial Time-Series Analysis
Here, (m) denotes the covariance at lag m for a block of data. There are di§erent ways
to compute (m), speciÖcally the scaling factor, so let us use the following convention:
assume the data is 0 mean; if not subtract the sample mean Örst. Then (with data indexed
1 n N):
(m) = 1
N
X
N
n=m+1
xnxnm
The (auto-)correlation coe¢ cients are:
(m) = (m) = (0)
For example, the MATLAB function autocorr(x) will compute and plot (m) for 0 m 20,
with signiÖcance bounds (in the graph, if a displayed value is below the bounds, then it can be
considered as 0). For our purposes, it is Öne if we use j (m)j > 0:2 as a test of ìsigniÖcanceî.
Our formula for (m) is technically biased, and can be interpreted as a Bartlett (i.e., similar
to triangular) windowed form of unbiased estimates. In any case, for our purposes, the lags
of interest are M, where M N so this distinction is not of interest for us, here.
For pairs of signals, we can deÖne:
xy (m) = 1
N
X
N
n=m+1
xnynm
for m 0, and similarly for m < 0. Then the cross-correlation coe¢ cient is:
xy (m) =
xy (m)
q
xx (0) yy (0)
and in MATLAB, crosscorr(x; y) will compute and graph this, by default for 20 m 20.
1. Heavy Tail Distributions
First we are going to explore some distributions. Generate N = 1e6 iid samples of
N (0; 1), Cauchy with = 1, and Studentsí t-distribution with = 5 and = 10
degrees of freedom. The variance of Studentsít-distribution is Önite for 3 (well,
technically > 2 but here we consider only integer values for ) and is given by:
2
For comparison, you should normalize your Studentsít-distribution data sets so they
have variance 1. Obviously donít normalize Cauchy since it has no variance! However,
= 0:544 is a reasonable “match” to N (0; 1) in the sense that both yield approximatel
the same P (jXj < 1) = 68% (which is the probability of being within one standard
deviation for the case of N (0; 1)):
In MATLAB, trnd generates the Studentsí t-distribution. You can generate Cauchy
= 1 as a special case of Studentsít by setting = 1; multiply the result by to
create Cauchy with general . Another approach is that if U is uniform over (0,1),
then X = tan (U) is Cauchy. Note in general Cauchy has pdf;
f (x) = =
2 + x
2
In any case, at this point you have N samples of ìcomparably scaledîGaussian, Cauchy
and Studentsí data. For each, compute the fraction of the time the absolute value
exceeds 4. Also, graph the Cauchy data: you will see something interesting (you will
not see in the others).
2. AR and ARMA
Let us start by synthesizing data using an ARMA(2; 2) model:
rt =
(1 + 0:2z
1
) (1 + 0:5z
1
)
(1 0:8z
1
) (1 0:7z
1
)
vt
where vt
is iid N (0; 1). Interpret the above as an operator notation for the appropriate
di§erence equation. Generate N = 250 samples. Consider the AR(p) model:
rt = vt + wp1rt1 + + wpprtp
The AR coe¢ cients as given in the notes would be ap0 = 1 and apk = wpk, 1 k p.
As usual we take ap0 = 1. We can interpret this as a linear regression model with model
error vt
. In this way, we can compute the AR coe¢ cients fapkg
p
k=1 using a least-squares
Öt. Recall app is the reáection coe¢ cient p, and E
jvt
j
2
is Pp, the power in the p
th
order prediction. The set p and Pp can be found using a Levinson-Durbin recursion
instead.
Take maximum order M = 10.
(a) Estimate the covariances (m) for 0 m M. Obtain a stem plot of the values
(m) = (0) for 0 m M. Those above say about 0:2 are signiÖcant. The
number above that level gives an indication that an AR model of that length
would be useful, higher order models would give marginal improvement. Note:
In MATLAB, autocorr does all this for you.
(b) Set up the corresponding (M + 1) (M + 1) Toeplitz matrix C, and compute its
eigenvalues. They are (hopefully) strictly positive real.
(c) Recall the Cholesky factorization of a pd matrix is C = LcholL
T
chol where Lchol is
lower triangular with positive entries on the diagonal. Closely related to this is
the so-called LDL decomposition C = LDLT where L is lower triangular with 1ís
on the diagonal, and D is a diagonal matrix with positive entries. If f`ig are the
diagonal elements of Lchol, then D = diag f`
2
i g. You should be able to
decomposition in your function library whether you use MATLAB, Python or
really any good scientiÖc computing software; if not and all you have is Cholesky,
you can derive D from Lchol as suggested, and then L = LcholD1=2
. If you canít
Önd either LDL or Cholesky decomposition, look harder. If itís not there, get a
better library! In any case, Önd the L and D for your covariance matrix C.
(d) Use the Levinson-Durbin recursion (you can use a function from your library, if
not then hand code it) to compute the reáection coe¢ cients and prediction error
powers Pm up to M. Also compute the FPEF coe¢ cients of all orders up to order
M. Pack these FPEFs into a triangular matrix of the form:
F =
2
6
6
6
6
6
4
1 0 0 0
a11 1 0 0
a22 a21 1 0
.
.
.
aMM aM1 1
3
7
7
7
7
7
5
(e) Compute F CFT
. Compare the diagonal elements of D with the Pmís:Compare
F with L
1
. In summary, the Levinson-Durbin recursion is tied to the LDL
decomposition (though it doesnít give you L and D directly, it gives what turns
out to be the more useful L
1 and D).
(f) Now use a LS Öt to compute the AR coe¢ cients of order M. Compare them with
the Mth order model you found via Levinson-Durbin.
(g) Comment on the reáection coe¢ cients you found. Consider the following general
issue: if a reáection coe¢ cient is close to 1 in magnitude, what does that tell you?
if it is close to 0, then what? Hint: Recall the formula that relates the Ppís to
the pís. Also, recall what holds for the ís if the data is exactly AR (M0) for
some order M0.
3. Repeat the above experiment with the modiÖcation that you have something close to
an ARIMA(2,1,2) model: r = Hv where:
H (z) =
1 0:99z
1
1 H0 (z)
where H0 (z) is the ARMA model in Problem 2. In other words, your rt will now exhibit
behavior that would suggest a unit-root nonstationarity. Plot the original rt (one
simulated block of data), and now the rt with this (close to) unit root nonstationarity.
After you generate the AR models and do the other analysis, now compute st = rtrt1,
the Örst di§erence, and repeat all for st
. The model should work better for st
. The
question is, can you detect the unit root nonstationarity early on? Go back to where
you graphed the partial correlation coe¢ cients, i.e., (m) = (0) for rt and comment.
4. Now repeat both Problems 2 and 3 using a Studentsít-distribution with = 5 for
vt (normalized so E (v
2
t
) = 1 still). The parts to repeat, and the issue of concern, is
the numerical stability and apparent Öt of the model. In theory, this should change
NOTHING! The theoretical values of (m) are IDENTICAL. There are various ways
of examining the changes, see if you can detect anything noticeable (or, in particular,
disturbing) happening. If not, good fo
5. Cointegration
Here we look for an example of cointegration, among other e§ects in time series. Let us
start with two underlying signals at
, bt
, with the Örst having a unit root nonstationarity
and the other stationary:
at
bt
= + G
at1
bt1
+
ut
vt
where:
=
0:3
0:5
G =
0:99 0
0:3 0:3
where ut
is ARMA(2; 2) driven by an input white noise that is iid N (0; 1), with innovations Ölter:
H =
(1 + 0:2z
1
) (1 + 0:3z
1
)
(1 0:2z
1
) (1 0:5z
1
)
and vt are iid samples with E (vt) = 0, E (v
2
t
) =
2
0
(a constant to be prescribed later),
and two possible distributions you will test: N (0; 2
0
), or Studentsít-distribution with
= 5 scaled to achieve the desired E (v
2
t
).
Next, let:
A =
3 4
2 3
and then:
rt =
r1t
r2t
= A
at
bt
(a) Compute the
0
vector and G0 matrix such that we can write:
r1t
r2t
=
0 + G
0
r1t1
r2t1
+
1t
2t
where 1t
; 2t are correlated disturbances, computed via t = A
ut
vt
.
(b) We want the e§ect of bt to be noticeable but not dominant, so do the following.
After you generate at
, take:
2
0 =
1
16
E
a
2
=
1
16
1
N
X
N
n=1
a
2
n
!
Then use that to generate bt
. Finally compute r1t
; r2t
. [In reality, you should do
this twice, using the two di§erent distributions for vt
. Do the rest of this exercise,
as well, for each case].
(c) Graph ut and vt
,. separately. Superimpose graphs of at
; bt
. Superimpose graphs
of r1t
; r2t
. Repeat this, say three times, to see “typical”
(d) We now want to detect the unit root stationarity in r1t
; r2t
, and determine if
there is indeed cointegration. Indeed, we can recover at
; bt via A1
rt
, which shows
that a linear combination of r1t
; r2t
is stationary! Determine the coe¢ cients of
cointegration (i.e., the linear combination of these two that yield a stationary
result).
(e) Compute and graph the auto-correlation coe¢ cients of r1t
; r2t
, separately, i.e.,
11 (m), 22 (m) for 0 m 20, and compute the cross-correlations 12 (m) for
20 m 20. In MATLAB, autocorr and crosscorr will do this for you.
(f) With sit = (1 z
1
) rit, i = 1; 2, repeat the above (i.e., look at the two autocorrelations and cross-correlation coe¢ cients between s1t
; s2t). You should see a
steeper decay which indicates both r1t
; r2t su§er from unit root nonstationarity.
(g) You will note the cross-correlation between r1t
; r2t decays slowly. To check this is
indicative of co-integration, try creating another signal that is INDEPENDENT
of r1t but also has a unit root nonstationarity. For example, let t be iid N (0; 1),
and ct = 0:99ct1 + t
. Checking the c
(m) coe¢ cients of c should of course
conÖrm it has the unit root nonstationarity. But now also look at the crosscorrelation coe¢ cients between ct and r1t
, and between ct and r2t
. You should see
signiÖcantly less cross-correlation between them.
6. ARCH/GARCH
Now we try ARCH/GARCH modeling. Take two years each of daily data for S&P500
and say two stocks that interest you, letís say some time within the last 20 years.
(You can get data o§ Yahoo Finance). SpeciÖcally, we are looking at adjusted closing
prices, and as a Örst step compute the daily log-returns for the period. Each year
will produce about 250 data points, but realize the exact number may vary.
First graph the returns and square returns. The returns should appear to be close to
constant mean (if not zero mean). Since we do not want to deal with the means here,
compute and subtract o§ the sample means, so from this point we will assume the
returns rt have 0 mean (constant). You may also try multiyear models, say 2007-2009,
which captures the 2008 Önancial crisis, or 2018-2020 (to date) which captures…. hmm,
that which will not be mentioned!
Compute and graph the auto-correlation coe¢ cients up to lag 20 for the returns AND
the square returns. If there are not signiÖcant values for the returns, that suggests our
mean model can be simple, i.e., constant mean. If the values for the square returns
however show some signiÖcance, then it is indicative a conditional heteroskedasticity
is in play. Ideally, for this simple experiment, your data should exhibit this (conditional heteroskedasticity and constant mean), and pick di§erent data sets if not. If we
assume the conditional expectation is 0 (letís say subtract o§ the sample mean Örst, if
necessary), our model is simply:
rt = t
t = tzt
2
t = GARCH (p; q)
; 2
where zt are iid with E (zt) = 0, E (z
2
t
) =
(a) First some simulation. For a GARCH(1; 1) model we can write:
2
t = ! + 2
t1 + 2
t1
Synthesize rt according to such a model, with zt N (0; 1), with ! = 0:5, = 0:6,
= 0:4. Superimpose graphs of rt and t
. You should see variability in the
volatility of rt
, that t seems to “track”. Now, repeat this with zt have a Students
t-distribution with = 5 degrees of freedom.
For each case, Öt a GARCH(1; 1) model, and an ARCH(2) model, assuming a
Gaussian distribution for zt (this is called QML- quasi-maximum likelihood, since
you are assuming a Gaussian likelihood function even when it is not the true one).
For the GARCH(1; 1) model, do you get back something close to your coe¢ cients?
For the ARCH(2) model, superimpose plots of rt and t as generated by the model,
and see if you have a good match for the volatility. What impact, if any, does the
Students t-distribution with = 5 degrees of freedom have?
(b) Your GARCH Ötting library should provide some information regarding the signiÖcance of the model coe¢ cients. Take your actual Önancial data, and Öt rt with
a GARCH(1,1) and ARCH(2) and see if either provides good Öts, and also plot rt
and t superimposed. If the volatility model is a good Öt, the graph of t should
appear to be like an ìenvelopeî for the return data.
