# GMMexamplesEcon583.ssc GMM examples for class slides # # authors: Eric Zivot and Jeff Wang # created: November 5, 2006 # revised: November 5, 2006 # # comments # 1. Notation follows Hayashi (2000) Econometrics, Princeton University Press # 2. consumption data in consump.ts is taken from Wooldridge (2002) # 3. CRSP data in pricing.ts is taken from Verbeek (2002) # module(finmetrics) # # 1. student t distribution # ref: Hamilton (1994) Time Series Analysis # chapter 15 # function to compute moment conditions for estimating # df parameter from Student-t distribution t.moments <- function(parm,data=NULL) { # parm = df parameter # data = [y^2, y^4] m1 = parm/(parm - 2) m2 = 3*parm*parm/((parm - 2)*(parm - 4)) t(t(data) - c(m1,m2)) } # simulate data from Student-t with 10 df set.seed(123) y = rt(250,df=10) y = y - mean(y) data = cbind(y^2,y^4) t.moments(10,data) colMeans(t.moments(10,data)) # show graphical analysis of t-distribution par(mfrow=c(2,2)) hist(y,probability=T) boxplot(y,outchar=T) plot(density(y),type="l",xlab="y(t)", ylab="density estimate") qqnorm(y) qqline(y) par(mfrow=c(1,1)) # estimate model by iterated efficient GMM start.vals = 15 names(start.vals) = c("theta") t.gmm.iter = GMM(start.vals,t.moments, method="iterative",max.steps=100, data=data) t.gmm.iter summary(t.gmm.iter) # sample moments evaluated at gmm estimate colMeans(t.moments(coef(t.gmm.iter),data)) # plot weight matrix pick.w = lower.tri(t.gmm.iter\$weight.matrix,diag=T) t.weights = t.gmm.iter\$weight.matrix[pick.w] barplot(t.weights, names=c("w11","w12","w22"), main="Efficient weights") # estimate model by CU efficient GMM start.vals = 15 names(start.vals) = c("theta") t.gmm.cu = GMM(start.vals,t.moments, method="simultaneous", data=data) t.gmm.cu summary(t.gmm.cu) # # 2. MA(1) model # ref: Harris (2000) # in Generalized Method of Moments, Matays (ed.), Cambridge # # function to compute four moments from MA(1) model # y(t) = mu + e(t) + psi*e(t-1) # e(t) ~ iid (0,sig2) ma1.moments <- function(parm,data=NULL) { # parm = (mu,psi,sig2)' # data = (y(t),y(t)^2,y(t)*y(t-1),y(t)*y(t-2)) is assumed to be a matrix m1 = parm[1] m2 = parm[1]^2 + parm[3]*(1 + parm[2]^2) m3 = parm[1]^2 + parm[3]*parm[2] m4 = parm[1]^2 t(t(data) - c(m1,m2,m3,m4)) } # simulate from MA(1) using arima.sim set.seed(123) ma1.sim = arima.sim(model=list(ma=-0.5),n=250) par(mfrow=c(3,1)) tsplot(ma1.sim,main="Simulated MA(1) Data") abline(h=0) tmp = acf(ma1.sim,plot=F) tmp2 = acf(ma1.sim,type="partial",plot=F) acf.plot(tmp,main="SACF") acf.plot(tmp2,main="SPACF") par(mfrow=c(1,1)) summaryStats(ma1.sim) # check moment function # data = (y(t),y(t)^2,y(t)*y(t-1),y(t)*y(t-2)) is assumed to be a matrix nobs = numRows(ma1.sim) ma1.data = cbind(ma1.sim[3:nobs],ma1.sim[3:nobs]^2, ma1.sim[3:nobs]*ma1.sim[2:(nobs-1)],ma1.sim[3:nobs]*ma1.sim[1:(nobs-2)]) start.vals = c(0,0.5,1) names(start.vals) = c("mu","psi","sig2") ma1.mom = ma1.moments(parm=start.vals,data=ma1.data) ma1.mom[1:5,] colMeans(ma1.mom) # check scaling, correlation and autocorrelations of moments at true parameters var(ma1.mom) cor(ma1.mom) tmp = acf(ma1.mom) start.vals = c(0,0.5,1) names(start.vals) = c("mu","psi","sig2") # estimate using truncated kernel with bandwith = 1 ma1.gmm.trunc = GMM(start.vals, ma1.moments, data=ma1.data,ts=T, var.hac.control=var.hac.control(bandwidth=1, window="truncated")) ma1.gmm.trunc summary(ma1.gmm.trunc) # sample moments evaluated at gmm estimate colMeans(t.moments(coef(t.gmm.iter),data)) # plot weight matrix pick.w = lower.tri(ma1.gmm.trunc\$weight.matrix,diag=T) ma1.weights = ma1.gmm.trunc\$weight.matrix[pick.w] barplot(ma1.weights, names=c("w11","w21","w31","w41","w22","w32","w42","w33","w43","w44"), main="Efficient weights") # MA(1) model fit to AR(1) data # y(t) = mu + u(t) # u(t) = phi*u(t-1) + e(t) # e(t) ~ iid (0,sig2) set.seed(123) ar1.sim = arima.sim(model=list(ar=0.5),n=250) par(mfrow=c(3,1)) tsplot(ar1.sim,main="Simulated AR(1) Data") abline(h=0) tmp = acf(ar1.sim,plot=F) tmp2 = acf(ar1.sim,type="partial",plot=F) acf.plot(tmp,main="SACF") acf.plot(tmp2,main="SPACF") par(mfrow=c(1,1)) summaryStats(ar1.sim) nobs = numRows(ar1.sim) ar1.data = cbind(ar1.sim[3:nobs],ar1.sim[3:nobs]^2, ar1.sim[3:nobs]*ar1.sim[2:(nobs-1)],ar1.sim[3:nobs]*ar1.sim[1:(nobs-2)]) start.vals = c(0,0.5,1) names(start.vals) = c("mu","psi","sig2") ar1.mom = ma1.moments(parm=start.vals,data=ar1.data) ar1.mom[1:5,] colMeans(ar1.mom) start.vals = c(0,0.5,1) names(start.vals) = c("mu","psi","sig2") # estimate using truncated kernel with bandwith = 1 ma1.gmm.ar1 = GMM(start.vals, ma1.moments, data=ar1.data,ts=T, var.hac.control=var.hac.control(bandwidth=1, window="truncated")) ma1.gmm.ar1 summary(ma1.gmm.ar1) # use one more moment # function to compute five moments from MA(1) model # y(t) = mu + e(t) + psi*e(t-1) # e(t) ~ iid (0,sig2) ma1.moments5 <- function(parm,data=NULL) { # parm = (mu,psi,sig2)' # data = (y(t),y(t)^2,y(t)*y(t-1),y(t)*y(t-2),y(t)*y(t-3)) m1 = parm[1] m2 = parm[1]^2 + parm[3]*(1 + parm[2]^2) m3 = parm[1]^2 + parm[3]*parm[2] m4 = parm[1]^2 m5 = parm[1]^2 t(t(data) - c(m1,m2,m3,m4,m5)) } nobs = numRows(ar1.sim) ar1.data = cbind(ar1.sim[4:nobs],ar1.sim[4:nobs]^2, ar1.sim[4:nobs]*ar1.sim[3:(nobs-1)],ar1.sim[4:nobs]*ar1.sim[2:(nobs-2)], ar1.sim[4:nobs]*ar1.sim[1:(nobs-3)]) start.vals = c(0,0.5,1) names(start.vals) = c("mu","psi","sig2") ar1.mom = ma1.moments5(parm=start.vals,data=ar1.data) ar1.mom[1:5,] colMeans(ar1.mom) start.vals = c(0,0.5,1) names(start.vals) = c("mu","psi","sig2") # estimate using truncated kernel with bandwith = 1 ma1.gmm.ar15 = GMM(start.vals, ma1.moments5, data=ar1.data,ts=T, var.hac.control=var.hac.control(bandwidth=1, window="truncated")) ma1.gmm.ar15 summary(ma1.gmm.ar15) # 3. Euler equation asset pricing model # ref: Hayashi (2000), Econometrics, Princeton University Press # Verbeek (2002), A Guide to Modern Econometrics, Wiley # Cochrane (2001), Asset Pricing, Princeton University Press # moment function for basic Euler equation asset pricing # model with power utility and one risky asset euler1.moments <- function(parm,data=NULL) { # parm = (beta,gamma) # data = (1+R(t+1),C(t+1)/C(t),1,C(t)/C(t-1), # C(t-1)/C(t-2),R(t),R(t-1),...) ncol = numCols(data) euler = data[,1]*parm[1]*(data[,2]^(-parm[2])) - 1 as.matrix(rep(euler,(ncol-2))*data[,3:ncol]) } # moment function for basic Euler equation asset pricing # model with power utility, J risky assets and 1 risk-free asset euler2.moments <- function(parm,data=NULL) { # parm = (beta,gamma) # data = (C(t+1)/C(t),1+Rf(t+1),R1(t+1)-Rf(t+1),..., # RJ(t+1)-Rf(t+1)) ncol = numCols(data) sdf = parm[1]*data[,1]^(-parm[2]) d1 = (sdf - 1)*data[,2] d2 = as.matrix(rep(sdf,(ncol-2))*data[,3:ncol]) cbind(d1,d2) } # Use pricing data from Verbeek for estimation pricing.ts@title pricing.ts@documentation colIds(pricing.ts) # plot pricing data plot(pricing.ts[,"CONS"], reference.grid=F, main="Consumption growth") tmp = seriesMerge(pricing.ts[,"CONS"],(1+pricing.ts[,-c(1,12)])) seriesPlot(tmp, one.plot=F, type="l") # create data series for examples pricing.mat = as.matrix(seriesData(pricing.ts)) ncol = numCols(pricing.mat) excessRet.mat = apply(pricing.mat[,2:(ncol-1)],2, function(x,y){x-y}, pricing.mat[,"RF"]) # data = (C(t+1)/C(t),1+Rf(t+1),R1(t+1)-Rf(t+1),..., # RJ(t+1)-Rf(t+1)) euler.data = cbind(pricing.mat[,"CONS"],1+pricing.mat[,"RF"], excessRet.mat) colIds(euler.data)[1] = "CONS" colIds(euler.data)[2] = "RF" # data = (1+R(t+1),C(t+1)/C(t),C(t)/C(t-1), # C(t-1)/C(t-2),R(t),R(t-1),...) n.obs = numRows(pricing.mat) euler.data1 = cbind(1+pricing.mat[3:n.obs,"R1"], pricing.mat[3:n.obs,"CONS"], rep(1,(n.obs-2)), pricing.mat[2:(n.obs-1),"CONS"], pricing.mat[1:(n.obs-2),"CONS"], 1+pricing.mat[2:(n.obs-1),"R1"], 1+pricing.mat[1:(n.obs-2),"R1"]) colIds(euler.data1) = c("1+R","CONS","CONST","CONS(-1)","CONS(-2)", "R(-1)","R(-2)") # # estimate Euler equation asset pricing model with # power utility and a single asset # # first few rows of euler1.moments start.vals = c(1,5) names(start.vals) = c("beta","alpha") euler1.mom = euler1.moments(start.vals,euler.data1) euler1.mom[1:5,] colMeans(euler1.mom) # gmm estimation # first do iterated GMM euler1.gmm.fit = GMM(start.vals, euler1.moments, method="iterative", data=euler.data1) euler1.gmm.fit summary(euler1.gmm.fit) # # estimate Euler equation asset pricing model with power utility, # J risky assets and a risk-free asset # # first few rows of euler2.moments start.vals = c(1,5) names(start.vals) = c("beta","alpha") euler.mom = euler2.moments(start.vals,euler.data) euler.mom[1:5,] colMeans(euler.mom) # gmm estimation # first do iterated GMM # note: results are very close to Verbeek pg. 146 start.vals = c(1,5) names(start.vals) = c("beta","alpha") euler.gmm.fit = GMM(start.vals, euler2.moments, method="iterative",data=euler.data) euler.gmm.fit summary(euler.gmm.fit) # next do inefficient GMM with W=Identity # results are very close to Verbeek pg. 146 # however, Verbeek uses non-standard J-stat euler.gmm.fit2 = GMM(start.vals, euler2.moments, method="iterative", max.steps=0, w=diag(11),w0.efficient=F, data=euler.data) euler.gmm.fit2 summary(euler.gmm.fit2) # compute average excess returns excessRet.hat = colMeans(excessRet.mat) # compute predicted excess returns from asset pricing model predRet <- function(parm,data) { # parm = (beta,gamma) # data = (C(t+1)/C(t),1+Rf(t+1),R1(t+1)-Rf(t+1),..., # RJ(t+1)-Rf(t+1)) ncol = numCols(data) sdf = parm[1]*data[,1]^(-parm[2]) tmp <- function(x,y) { -var(x,y)/mean(y) } ans = apply(data[,3:ncol],2,FUN=tmp,sdf) ans } predRet.hat = predRet(coef(euler.gmm.fit2),euler.data) # plot avgerage excess returns against predicted returns plot(predRet.hat,excessRet.hat, ylab="Mean excess return", xlab="Predicted mean excess return", xlim=c(0,0.01),ylim=c(0,0.01)) abline(a=0,b=1) # annualize monthly data to match Verbeek pg. 147 predRet.hat.a = (1+predRet.hat)^12 - 1 excessRet.hat.a = (1+excessRet.hat)^12 - 1 plot(predRet.hat.a,excessRet.hat.a,ylab="Mean excess return", xlab="Predicted mean excess return",xlim=c(0,0.15),ylim=c(0,0.15)) abline(a=0,b=1) # # 4. Stochastic Volatility model # ref: Andersen and Sorensen (1996) # ``GMM Estimation of a # Stochastic Volatility Model: A Monte Carlo Study'', # Journal of Business and Economic Statistics # # use simulator function for AS SV model # from EMM chapter # # y(t) = exp(w(t)/2)*z(t) # w(t) = alpha + beta*w(t-1) + sigma.u*u(t) # -1 < beta < 1, sigma.u > 0 # simulator function # note: need transformation to impose stationarity on simulations # set beta = exp(beta.t)/(1+exp(beta.t)) # beta.t = log(beta/(1-beta)) sv.as.gensim <- function(rho,n.sim,n.var=1,n.burn,aux=sv.as.aux) { # rho = (alpha, beta.t, sigma.u) # aux is a list with components # z = vector of N(0,1) values of length (n.sim + n.burn) # u = vector N(0,1) values of length (n.sim + n.burn) nz = n.burn + n.sim beta = exp(rho[2])/(1+exp(rho[2])) mu = rho[1]/(1-beta) w = mu + arima.sim(model = list(ar=beta), innov = aux\$u[(n.burn+1):nz]*rho[3], start.innov = aux\$u[1:n.burn]*rho[3]) y = exp(w/2)*aux\$z[(n.burn+1):nz] as.numeric(y) } # simulate data from AS Monte Carlo design 1 # omega =-0.736, beta=0.90, sigma.u=0.363 # beta.t = log(beta/(1-beta)) n.sim = 1000 n.burn = 100 nz = n.sim+n.burn set.seed(456) sv.as.aux = list(z=rnorm(nz),u=rnorm(nz)) rho.as = c(-0.736,2.197225,0.363) names(rho.as) = c("alpha","beta.t","sigma.u") sv.as.sim = sv.as.gensim(rho=rho.as,n.sim=n.sim, n.burn=n.burn,aux=sv.as.aux) # plot simulated data tsplot(sv.as.sim) abline(h=0) summaryStats(sv.as.sim) par(mfrow=c(2,1)) tmp = acf(sv.as.sim) tmp = acf(abs(sv.as.sim)) par(mfrow=c(1,1)) par(mfrow=c(2,2)) hist(sv.as.sim,probability=T) boxplot(sv.as.sim,outchar=T) plot(density(sv.as.sim),type="l",xlab="y(t)", ylab="density estimate") qqnorm(sv.as.sim) qqline(sv.as.sim) par(mfrow=c(1,1)) ## define moment equations sv.moments = function(parm, data=NULL) { omega = parm[1] beta = parm[2] sigu = parm[3] mu = omega/(1-beta) sig2 = (sigu*sigu)/(1-beta*beta) # E.sigma = c(sqrt(2/pi) * exp(mu/2 + sig2/8), exp(mu + sig2/2), 2 * sqrt(2/pi) * exp(3*mu/2 + 9*sig2/8), 3 * exp(2*mu + 2*sig2)) E.sigma.c = c(2/pi * exp(2*(mu/2 + sig2/8) + beta^(1:10) * sig2/4), exp(2*(mu + sig2/2) + 4 * beta^(1:10) * sig2/4)) # t(t(data) - c(E.sigma, E.sigma.c)) } # create data to be passed to sv.moments function sv.pow = cbind(abs(sv.as.sim),sv.as.sim^2,abs(sv.as.sim)^3, sv.as.sim^4) sv.pow = sv.pow[-(1:10),] sv.cm = tslag(sv.as.sim, 1:10, trim=T) * as.vector(sv.as.sim[-(1:10)]) sv.data = cbind(sv.pow, abs(sv.cm), sv.cm^2) colIds(sv.data) = NULL # iterative GMM estimation of the SV model start.vals = c(0,0.5,0.5) names(start.vals) = c("omega","beta","sigu") sv.fit.1 = GMM(start.vals, sv.moments, method="iterative", ts=T, data=sv.data) summary(sv.fit.1) # CU GMM estimation sv.fit.2 = GMM(c(0.4, 0.5, 0.5), sv.moments, method="simul", ts=T, data=sv.data) sv.fit.3 = GMM(c(0.4, 0.5, 0.5), sv.moments, method="simul", ts=T, data=sv.data, var.hac.control=var.hac.control(automatic="nw")) # fit SV model to actual sp500 data plot(sp500.ts, reference.grid=F) abline(h=0) summaryStats(sp500.ts) SP500.ts = sp500.ts SP500.ts@data = as.matrix(seriesData(sp500.ts)) SP500.pow = seriesMerge(abs(SP500.ts), SP500.ts^2, abs(SP500.ts)^3, SP500.ts^4) SP500.pow = SP500.pow[-(1:10),]@data SP500.cm = tslag(SP500.ts, 1:10, trim=T)@data * as.vector(SP500.ts[-(1:10)]@data) SP500.data = cbind(SP500.pow, abs(SP500.cm), SP500.cm^2) colIds(SP500.data) = NULL # iterative GMM estimation of the SV model start.vals = c(0,0.5,0.5) names(start.vals) = c("omega","beta","sigu") sv.fit.SP500 = GMM(start.vals, sv.moments, method="iterative", ts=T, data=SP500.data) summary(sv.fit.SP500,print.moments=F) # # 5. Interest rate diffusion model # ref: Chan, Karolyi, Longstaff and Sanders (1992) # Journal of Finance # James and Webber (2000). Interest Rate Models # # function to compute CKLS moments # note: dt defines the discretization step. For example, # if r(t) is annualized rate sampled monthly, set dt = 1/12 # if r(t) is monthly rate sampled monthly, set dt = 1 ckls.moments <- function(parm,data=NULL,dt=1/12) { # parm = (alpha,beta,sigma,gamma)' # data = [r(t+dt)-r(t),r(t)] # dt = discretization step e.hat = as.vector(data[,1] - (parm[1] + parm[2]*data[,2])*dt) m2 = e.hat*as.vector(data[,2]) m3 = e.hat^2 - dt*parm[3]*parm[3]*(as.vector(data[,2])^(2*parm[4])) m4 = m3*data[,2] cbind(e.hat,m2,m3,m4) } # # use Mike Cliff's data # td = timeSeq(from="6/30/1964",to="11/30/1989",by="months", format="%b %Y") ckls.ts = timeSeries(pos=td,data=ckls.df\$yield)/100 colIds(ckls.ts) = "RF" ckls.ts@title = "Annualized 30-day T-bill Rate" ckls.ts@documentation = "source: http://www.mgmt.purdue.edu/faculty/mcliff/progs.html" # compute data for GMM estimation data.ckls.ts = seriesMerge(diff(ckls.ts),tslag(ckls.ts)) colIds(data.ckls.ts)[1] = "RF.diff" data.ckls = as.matrix(seriesData(data.ckls.ts)) # summarize data par(mfrow=c(2,1)) plot(ckls.ts,reference.grid=F) plot(diff(ckls.ts),reference.grid=F) abline(h=0) par(mfrow=c(1,1)) summaryStats(data.ckls.ts) # note: set dt = 1/12 because data are annualized rate sampled monthly start.vals = c(0.04,-0.6,sqrt(1.7),1.5) ckls.mom = ckls.moments(parm=start.vals,data=data.ckls,dt=1/12) ckls.mom[1:5,] colMeans(ckls.mom) # GMM estimation of just identified ckls model start.vals = c(0.06,-0.5,1,1) names(start.vals) = c("alpha","beta","sigma","gamma") gmm.ckls = GMM(start.vals,ckls.moments,ts=T, data=data.ckls,dt=1/12) summary(gmm.ckls) # compute long-run mean theta.hat = coef(gmm.ckls) theta.hat["alpha"]/-theta.hat["beta"] -theta.hat["beta"] # Estimate restricted model cir.moments <- function(parm,data=NULL,dt=1/12) { # parm = (alpha,beta,sigma)' # data = [r(t+dt)-r(t),r(t)] # dt = discretization step e.hat = as.vector(data[,1] - (parm[1] + parm[2]*data[,2])*dt) m2 = e.hat*as.vector(data[,2]) m3 = e.hat^2 - dt*parm[3]*parm[3]*(as.vector(data[,2])) m4 = m3*data[,2] cbind(e.hat,m2,m3,m4) } # iterated GMM estimation of overidentified CIR model start.vals = c(0.06,-0.5,1) names(start.vals) = c("alpha","beta","sigma") gmm.cir = GMM(start.vals,cir.moments,ts=T, data=data.ckls,dt=1/12) summary(gmm.cir) # compute long-run mean theta.hat = coef(gmm.cir) theta.hat["alpha"]/-theta.hat["beta"] -theta.hat["beta"]