# GMMexamples.ssc examples for GMM documentation # # authors: Eric Zivot and Jeff Wang # created: November 1, 2003 # revised: November 19, 2004 # revised: March 23, 2005 # # 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) # ########################################################################## # Estimation of S ########################################################################## args(var.hac) # function to compute kernel weights kernel.weights <- function(n,lags) { bn.trunc = trunc(4*(n/100)^(1/5)) bn.bartlett= trunc(4*(n/100)^(1/4)) bn.parzen = trunc(4*(n/100)^(4/25)) bn.qs = trunc(4*(n/100)^(4/25)) ai.trunc = seq(1:lags)/(1+bn.trunc) ai.bartlett= seq(1:lags)/(1+bn.bartlett) ai.parzen = seq(1:lags)/(1+bn.parzen) di.qs = seq(1:lags)/(bn.qs) mi.qs = 6*pi*di.qs/5 wt.trunc = ifelse((ai.trunc < 1), 1, 0) wt.bartlett = ifelse((ai.bartlett<= 1), 1-ai.bartlett, 0) wt.parzen = rep(0,lags) wt.parzen = ifelse((ai.parzen <= 0.5),1 - 6*ai.parzen^2 + 6*ai.parzen^3, 2*(1 - ai.parzen)^3) wt.parzen[(ai.parzen > 1)] = 0 wt.qs = (25/(12*di.qs)^2)*(sin(mi.qs)/mi.qs - cos(mi.qs)) ans = cbind(wt.trunc,wt.bartlett,wt.parzen,wt.qs) ans } # plot kernel weights k.wts = kernel.weights(n=100,lags=10) tsplot(k.wts, ylab = "weights", xlab="lags", lty=1:4,lwd=2) legend(6,1,legend=c("Truncated","Bartlett","Parzen","Quadratic Spectral"), lty = 1:4,lwd=2) # compute long-run variance for bivariate time series dspot = diff(lexrates.dat[,"USCNS"]) fp = lexrates.dat[,"USCNF"]-lexrates.dat[,"USCNS"] uscn.ts = seriesMerge(dspot,fp)*100 colIds(uscn.ts) = c("dspot","fp") uscn.ts@title = "US/CN Exchange Rate Data" par(mfrow=c(2,1)) plot(uscn.ts[,"dspot"],main="1st difference of US/CN spot exchange rate") plot(uscn.ts[,"fp"],main="US/CN interest rate differential") par(mfrow=c(1,1)) # compute short run cov var.hat = var(uscn.ts) var.hat # compute long-run variance S.hat.b = var.hac(uscn.ts, window="bartlett", bandwidth=4, demean=T) S.hat.b S.hat.default = var.hac(uscn.ts) S.hat.nw = var.hac(uscn.ts, window="bartlett", automatic="nw", prewhiten = F, w = c(1,1), demean = T) S.hat.am = var.hac(uscn.ts, window="qs", automatic="andrews", prewhiten = F, w = c(1,1), demean = T) S.hat.qs = var.hac(uscn.ts, window="qs", prewhiten = T, demean = T) S.hat.b S.hat.default S.hat.nw S.hat.am S.hat.qs # GMM arguments args(GMM) ########################################################################## # linear GMM ########################################################################## # # 1. OLS as GMM # # extract data from timeSeries excessReturns.df = seriesData(excessReturns.ts) # ols estimation of CAPM regression ols.fit = OLS(MSFT~SP500,data=excessReturns.ts) ols.fit # OLS moment conditions for GMM estimation ols.moments = function(parm,y=NULL,x=NULL) { x = as.matrix(x) x*as.vector(y - x%*%parm) } # test moment condition function ols.moments(c(1,1),y=excessReturns.df[,"MSFT"], x=cbind(1, excessReturns.df[,"SP500"])) # ols estimation using GMM # compare to OLS with White HC standard errors start.vals = c(0,1) names(start.vals) = c("alpha","beta") gmm.fit = GMM(start.vals, ols.moments, max.steps=1, y=excessReturns.df[,"MSFT"], x=cbind(1, excessReturns.df[,"SP500"])) class(gmm.fit) names(gmm.fit) gmm.fit summary(gmm.fit) summary(ols.fit,correction="white") # ols estmation using GMM with HAC covariance matrix # this gives the same result as summary(ols.fit, correction="nw") gmm.fit2 = GMM(c(0,1), ols.moments, max.steps=1, y=excessReturns.df[,"MSFT"], x=cbind(1, excessReturns.df[,"SP500"]), ts=T, df.correction = F, var.hac.control=var.hac.control(window="bartlett", bandwidth=floor(4 * (nrow(excessReturns.df)/100)^(2/9)))) summary(gmm.fit2) summary(ols.fit,correction="nw") # use df.correction=F to match Eviews gmm.fit2 = GMM(c(0,1), ols.moments, max.steps=1, y=excessReturns.df[,"MSFT"], x=cbind(1, excessReturns.df[,"SP500"]), ts=T, df.correction = F, var.hac.control=var.hac.control(window="bartlett", bandwidth=floor(4 * (nrow(excessReturns.df)/100)^(2/9)))) summary(gmm.fit2) # check with pre-whitened residuals # doesn't work gmm.fit2 = GMM(c(0,1), ols.moments, max.steps=1, y=excessReturns.df[,"MSFT"], x=cbind(1, excessReturns.df[,"SP500"]), ts=T, df.correction = F, var.hac.control=var.hac.control(window="bartlett", prewhiten=T, bandwidth=floor(4 * (nrow(excessReturns.df)/100)^(2/9)))) summary(gmm.fit2) # check with automatic bandwidth procedures gmm.fit2 = GMM(c(0,1), ols.moments, max.steps=1, y=excessReturns.df[,"MSFT"], x=cbind(1, excessReturns.df[,"SP500"]), ts=T, df.correction = F, var.hac.control=var.hac.control(window="bartlett", prewhiten=F, automatic="nw",w=c(0,1))) summary(gmm.fit2) # # 2. IV as GMM # # create moment function for IV estimation by GMM # y(t) = z(t)*delta + e(t), E[z(t)e(t)] /= 0 # E[x(t)e(t)] = 0 iv.moments = function(parm,y,X,Z) { X = as.matrix(X) Z = as.matrix(Z) X*as.vector(y - Z%*%parm) } # Campbell and Mankiw (1990) PIH example # using consumption data from Wooldridge (2002) # consump.ts@title = c("Consumption data from Wooldridge") # consump.ts@documentation = c("GY = log growth rate of real income", # "GC = log growth rate of real consumption", # "R3 = real 3-month T-bill rate", # "source: Wooldridge (2002), Introduction to", # "Econometrics, 2nd Edition", # "South-Western Thompson") colIds(consump.ts) plot(consump.ts,reference.grid=F, plot.args=list(lty=1:3)) legend(0,-0.01,legend=colIds(consump.ts),lty=1:3) # create data frame for analysis - add constant nobs = numRows(consump.ts) consump = seriesData(consump.ts) consump\$const = rep(1,nobs) # check moment condition function iv.mom = iv.moments(c(1,1,1),y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) iv.mom[1:5,] # one step GMM with a fixed inefficient weight matrix # W = I start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.1step = GMM(start.vals,iv.moments, method="iterative", max.steps=0, w=diag(4), w0.efficient=F, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) summary(gmm.fit.1step) # two step efficient estimator using identity matrix # as initial weight matrix # note: result won't match Eviews because Eviews uses # X'X as initial weight matrix and df correction start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.2step = GMM(start.vals,iv.moments, method="iterative", max.steps=1, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) summary(gmm.fit.2step) # two step efficient estimator using Sxx^-1 as weight matrix # matches Eviews coefficient output with max iterations = 1 w.tsls = crossprod(consump[1:(nobs-1),])/nobs start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.2step.eviews = GMM(start.vals,iv.moments, method="iterative",max.steps=1, w=w.tsls,w0.efficient = T,df.correction=F, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) # note: J-statistic and avar are not computed correctly because weight # matrix does not contain TSLS estimate of error variance summary(gmm.fit.2step.eviews) # iterated efficient estimator start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.iter = GMM(start.vals,iv.moments, method="iterative", max.steps=100, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) summary(gmm.fit.iter) # iterative estimator with df.correction=F to match Eviews output start.vals = rep(0.1,3) gmm.fit.eviews = GMM(start.vals,iv.moments, method="iterative",df.correction=F, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) summary(gmm.fit.eviews) # try simultaneous estimation - good starting values are important! start.vals = gmm.fit.iter\$parameters gmm.fit.cu = GMM(start.vals,iv.moments, method="simultaneous", y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) summary(gmm.fit.cu) # do gmm estimation with fixed weight matrix equal to the # efficient weight matrix from the iterated estimator # note: you need to invert the weight matrix to get the correct result w.fixed = solve(gmm.fit.iter\$weight.matrix) # use Jeff's trick to inforce symmetry w.fixed = (w.fixed + t(w.fixed))/2 start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.fixed = GMM(start.vals,iv.moments, method="iterative", max.steps=0, w=w.fixed,w0.efficient=T, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) summary(gmm.fit.fixed) summary(gmm.fit.iter) # # 3. do tsls as GMM estimator # # efficient weight matrix is Sxx^-1 w.tsls = crossprod(consump[1:(nobs-1),])/nobs start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.tsls = GMM(start.vals,iv.moments, method="iterative",max.steps=0, w=w.tsls,w0.efficient = T, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) # J-statistic and avar are not computed correctly because weight # matrix does not contain TSLS estimate of error variance # Why is the J-statistic # so close to zero? This because the error variance is not in the # weight matrix gmm.fit.tsls # compute TSLS estimate of error variance y = as.vector(consump[2:nobs,1]) X=as.matrix(consump[1:(nobs-1),]) Z=as.matrix(consump[2:nobs,2:4]) d.hat = coef(gmm.fit.tsls) e.hat = y - Z%*%d.hat df = nrow(Z) - ncol(Z) s2 = as.numeric(crossprod(e.hat)/df) # compute correct efficient weight matrix for tsls # that contains error variance term w.tsls2 = crossprod(X)*s2/nobs start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.tsls2 = GMM(start.vals,iv.moments, method="iterative",max.steps=0,df.correction=F, w=w.tsls2,w0.efficient=T, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) summary(gmm.fit.tsls2) # GMM converges faster with X'X used as initial weight matrix # note: why is 1-step objective equal to zero???? w.tsls = crossprod(consump[1:(nobs-1),])/nobs start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.iter2 = GMM(start.vals,iv.moments, method="iterative",max.steps=100, w=w.tsls,w0.efficient=F, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) # # 4. hypothesis testing # # test PIH: d1=d2=0 using Wald statistic Rmat = matrix(c(1,0,0,1,0,0),2,3) rvec = c(0,0) dhat = coef(gmm.fit.iter) avarRbhat = Rmat%*%gmm.fit.iter\$vcov%*%t(Rmat) Rmr = Rmat%*%dhat - rvec wald.stat = as.numeric(t(Rmr)%*%solve(avarRbhat)%*%Rmr) wald.stat 1 - pchisq(wald.stat,2) # test PIH: d1=d2=0 using GMM LR statistic # first compute restricted GMM estimator s.ur = solve(gmm.fit.iter\$weight.matrix) s.ur = (s.ur + t(s.ur))/2 start.vals = 0.1 names(start.vals) = c("const") gmm.fit.r = GMM(start.vals,iv.moments, method="iterative", max.steps=0, w=s.ur, w0.efficient=T, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=consump[2:nobs,4]) summary(gmm.fit.r) # statistic should be equivalent to wald statistic gmm.lr = gmm.fit.r\$objective - gmm.fit.iter\$objective gmm.lr 1 - pchisq(gmm.lr,2) # compare Wald and GMM-LR for nonlinear hypothesis # d1/d2 = 1 # restricted GMM fit # restricted model is # dc = d1 + d2*(dy + r) + error s.ur = solve(gmm.fit.iter\$weight.matrix) s.ur = (s.ur + t(s.ur))/2 start.vals = c(0.1,0.1) names(start.vals) = c("const","GY+R3") gmm.fit.r = GMM(start.vals,iv.moments, method="iterative", max.steps=0, w=s.ur, w0.efficient=T, y=consump[2:nobs,1], X=consump[1:(nobs-1),], Z=cbind(consump[2:nobs,4], (consump[2:nobs,2]+consump[2:nobs,3]))) summary(gmm.fit.r) gmm.lr = gmm.fit.r\$objective - gmm.fit.iter\$objective gmm.lr 1 - pchisq(gmm.lr,2) # # 5. Compute C stat for endogeneity of r(t) # # Full model assuming r(t) is exog # iterated efficient estimator start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.full = GMM(start.vals,iv.moments, method="iterative", max.steps=100, y=consump[2:nobs,"GC"], X=cbind(consump[1:(nobs-1),],consump[2:nobs,"R3"]), Z=consump[2:nobs,2:4]) summary(gmm.fit.full) # extract k1 x k1 component of efficient weight matrix w11.full = solve(gmm.fit.full\$weight.matrix[1:4,1:4]) w11.full = (w11.full + t(w11.full))/2 # reduced model assuming r(t) is endogenous estimated using # inverse of S11.full as efficient weight matrix start.vals = rep(0.1,3) names(start.vals) = c("GY","R3","const") gmm.fit.11 = GMM(start.vals,iv.moments, method="iterative", max.steps=0, w=w11.full, w0.efficient=T, y=consump[2:nobs,"GC"], X=consump[1:(nobs-1),], Z=consump[2:nobs,2:4]) summary(gmm.fit.11) # compute c statistic C.stat = gmm.fit.full\$objective - gmm.fit.11\$objective C.stat 1 - pchisq(C.stat,1) # # instrument relevance tests # # first stage regressions for pih model firstStage.fit = OLS(cbind(GY,R3)~tslag(GC)+tslag(GY)+tslag(R3), data=consump) class(firstStage.fit) # note: bug in summary.mOLS summary(firstStage.fit) ########################################################################## # nonlinear GMM ########################################################################## # # 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) # 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) # 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) # 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) # 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) ## 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) 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"]