MS&E349 Homework 2 Solution

Please submit one homework assignment per group. The homework solution should be submitted online to Canvas. Please indicate clearly the names of all group members. I prefer that solutions are typed in Latex, but it is also fine to submit scanned copies of handwritten solutions. Include the commented code in an appendix section. Please also submit the executable and commented code.
Theoretical Questions
Question 1 Method of Moments
Consider a random variable X that has an Poisson distribution with parameter λ, i.e. the probability of X is
for k ≥0 and P(X = k) = 0 for k < 0
Assume you observe a sequence of i.i.d. observations of random variables with Poisson distribution
{Xt : 1 ≤ t ≤ T}. The first four moments are given by
E[X] = λ E[X2] = λ2 + λ E[X3] = λ3 + 3λ2 + λ E[X4] = λ4 + 6λ3 + 7λ2 + λ.
1. (Bonus question: not required and no additional points) Calculate the first four moments of X.
(Hint: you can use and k4 = k(k −1)(k −2)(k −3)+6k(k −1)(k −2)+7k(k −1)+k to calculate the moments.)
2. Derive the method of moment estimator for λ using the first moment.
3. Derive the method of moment estimator for λ using the second moment.
4. Now consider the GMM estimator for λ using the first and second moment. Calculate the asymptotic distribution of λˆGMM. What is the optimal weighting matrix for the GMM estimator? How does the optimal GMM estimator compare to the estimators based only on one moment.
.........
Question 2 GARCH
Consider the GJR-GARCH model

Assume that the process is covariance stationary.
1. What is E[rt+1]?
2. What is Et[rt+1]?
3. What is V ar(rt+1)?
4. What is V art(rt+1)?
5. What is V art(rt+2)?
.........
Empirical Questions
The data that you need for this exercise can be found in the CRSP data base of WRDS. You have free access to CRSP with your Stanford ID. Please use the Stanford library link:
http://libguides.stanford.edu/databases/crsp.
Please download the following data:
Question 3 Volatility
In this question we want to estimate the volatility of the S&P 500 index. Consider a GARCH(1,1) model of the form

1. Simulate a time series of daily log returns for Tday = 13,200 days of a GARCH(1,1) model with µ = 0.01,ω = 0.0073,α = 0.93 and β = 0.06. Aggregate the daily log returns to Tmonth = 600 monthly returns by adding up 22 daily returns. This will be your simulated daily and monthly data set. The S&P 500 daily and monthly returns are the second data set.
2. Use the monthly simulated and return data of the S&P 500 index to estimate a GARCH(1,1) for monthly volatility. What features of the model and/or data are noteworthy? Report your estimates. You will need to use a recursive formulation in your code. Assume that all observations needed for starting the recursion are equal to their unconditional expectation.
3. Use the daily data (simulated and S&P 500 returns) to calculate a monthly time-series of realized volatilities (based on quadratic variations). (Assume that a month has 22 trading days.)
4. Compare your realized volatilities with your estimates based on the GARCH(1,1) model (for simulated and S&P 500 returns). Plot both time-series into the same graph.
5. Estimate an AR(1) model for realized volatility (for simulated and S&P 500 returns). What features of the model and/or data are noteworthy?
6. Use the estimated monthly realized volatilities and the estimated AR(1) model to produce one-month-ahead forecasts of the volatilities. Plot realized volatility and forecasted volatility.
7. Compute the Mean-Square-Error (MSE) of the AR(1) forecasts, defined as E[(σt+1−σˆt+1|t)2], where σt+1 is the realized value of the volatility at t + 1 and ˆσt+1|t is its forecast given information at time t.
8. Use the estimated monthly volatilities and the estimated GARCH(1,1) model to produce onemonth-ahead forecasts of the volatilities. Plot volatility as estimated by the GARCH(1,1) and forecasted volatility based on the GARCH(1,1).
9. Compute the Mean-Square-Error (MSE) of the GARCH(1,1) forecasts, defined as E[(σt+1 − σˆt+1|t)2], where σt+1 is the realized volatility based on the GARCH(1,1) model at t + 1 and σˆt+1|t is its forecast given information at time t based on the GARCH(1,1).
.........
Question 4 Parameter Estimation of SDEs
In this problem we estimate the parameters of the Cox- Ingersoll-Ross (CIR) process, which is often used to model interest rates and is defined as:
√
drt = κ(θ − rt)dt + σ rtdWt (1)
The model can also be used to model the volatility in stochastic variance models, since it always remain positive. The estimation will be applied to two data sets: (1) simulated data and (2) interest rate data.
Use an Euler approximation to obtain the following discrete time model of the CIR:
(2)
where ∆k = tk+1 − tk and is normally distributed with

and Ftk denotes the information set at time tk.
Simulate 100,000 steps with the Euler discretization and take every 250 observation of the process to obtain a simulated CIR model with 400 observations. Use the parameters κ = 0.234,θ = 0.081 and σ2 = 0.0073. Use the same simulated time series with 400 observations for all the estimations. The second time series is based on the monthly treasury yields. Report all your estimation results for both time series. The following questions will ask you to estimate the parameters
with different estimation approaches and to report the asymptotic covariance ma-
√ d
trix Σ from T(γˆ − γ) → N(0,Σ). Please report also the standard errors for each estimate ˆγi, i.e. . (Hint: You can calculate derivatives numerically.)
1. GMM: Estimate the parameters in the CIR model and the covariance matrices of the estimators by using GMM. For this purpose we approximate the stochastic differential equation with an Euler scheme (higher order methods would work even better). The advantage of using the discretized process is that we can estimate complex models without much problem. Use the following four moment restrictions:

In order to obtain an efficient GMM estimation, you can first estimate the parameters with the identity matrix as weighting matrix. Use these parameters to calculate the optimal weight matrix, chosen as the inverse of the covariance matrix for the moment restrictions. Estimate the new parameters with the new weight matrix.
2. QML: Estimate the parameters and their covariance matrices in the models above by using approximate Quasi Maximum Likelihood (by discretizing the model and using the likelihood generated by the approximate model).
3. Transformed approximated likelihood: Use a transformation of the data for the CIR-
√ model and the Shoji and Ozaki approximated likelihood. Using the transformation yt = 2 rt we obtain a state independent diffusion term. The Ito formula yields

We can now use Shoji and Ozaki’s approach to approximate the likelihood with a Gaussian likelihood:

4. Exact likelihood for the CIR model: The CIR process has a known transition density.
If we let ∆ , and Y = 2cX, then Ytk|Ytk−1 is distributed as a noncentral chi-squared with 2κθ/σ2 degrees of freedom and noncentrality parameter Ytk−1 exp(κ∆k). Using this we can write up the density as

where u = cxtk−1 exp(−κ∆k),v = cxtk,q = 2κθ/σ2 − 1 and Iq(z) is the modified Bessel function of the first kind of order q (iv(v,z) in Python or besseli(q,z) in Matlab).
Estimate the parameters and their covariance matrix for the Cox-Ingersoll-Ross model by using the exact likelihood.
5. Simulated Maximum Likelihood (Bonus question: Not required and no additional points)
Sometimes it is possible to use Maximum Likelihood on a simulated likelihood in stochastic differential equations. Assume that we have observations yn,n = 1,...,N from some model.
Also assume that the sampling interval is ∆. We now want to find the likelihood
N
L(θ) = Yf(yi|yi−1,θ).
i=1
If this is not available in closed form, one might approximate it using simulations. First we discretize the dynamics of the stochastic differential equation using a scheme that gives a simple transition probability, preferably Euler. Then, divide the interval ∆ into M subintervals of length δ = ∆/M. The idea is the to simulate K trajectories on a grid of size δ up to subinterval M − 1 starting in yn. Do this for every time n, resulting in K samples at every n. Under the Euler discretization, the transition density from M −1 to M is Gaussian with mean µn,k and standard deviation σn,k , both given by the model in question. Then we approximate the likelihood by the following:
!
where φ(y,µ,σ) is the density of the Normal distribution with mean µ and standard deviation σ. The log likelihood is given by

We then use numerical optimization techniques to maximize the loglikelihood. In this exercise we keep the number of intermediate steps fairly small, say M = 2 or 3. If M is too big, the variance of the sample will become large and make the estimation difficult. The remedy for this is to use an importance sampler.
It is usually a good idea to use the same sequence of random numbers each time the likelihood is evaluated. This is referred to as Common Random Numbers and is a way to avoid the Monte Carlo error in the minimization. Practically this is done by drawing a large number of random numbers, in this case N · (M − 1) · K and use these as input to the likelihood.
Estimate the parameters and their covariance matrices in the models above by using Simulated Maximum Likelihood. Choose an appropriate number of samples K based on your estimation results.
6. Compare your results for the different estimators.
.........