Chapter 11. Method ofMoments#
Empirical Results in Chapter 11#
dat <- read.table("cps09mar.txt")
experience <- dat[,1]-dat[,4]-6
mwm <- (dat[,11]==1)&(dat[,12]<=2)&(dat[,2]==0)&(experience==20)
dat1 <- dat[mwm,]
wage <- as.matrix(dat1[,5]/(dat1[,6]*dat1[,7]))
lwage <- log(wage)
education <- as.matrix(dat1[,4])
n <- length(wage)
mwage <- mean(wage)
mlwage <- mean(lwage)
meducation <- mean(education)
vwage <- var(wage)
vlwage <- var(lwage)
veducation <- var(education)
x <- cbind(wage,education)
c <- cov(x)
swage <- sd(wage)
slwage <- sd(lwage)
seducation <- sd(education)
sewage <- swage/sqrt(n)
seeducation <- seducation/sqrt(n)
wwage <- wage/100
wagem <- (wwage%*%matrix(1,1,4))^(matrix(1,n,1)%*%(1:4))
mwagem <- colMeans(wagem)
semwagem <- sqrt(colMeans((wagem - matrix(1,n,1)%*%mwagem)^2)/(n-1))
cwage <- wwage - mean(wwage)
cwagem <- (cwage%*%matrix(1,1,5))^(matrix(1,n,1)%*%(1:5))
cmwagem <- colMeans(cwagem)
cmwagem[1] <- mwagem[1]
v <- cov(cwagem)
secmwage <- matrix(0,1,4)
secmwage[1] <- var(wwage)
secmwage[2] <- var(cwage^2)
for (m in 3:4) {
secmwage[m] <- v[m,m] + ((m*cmwagem[m-1])^2)*var(wwage) -2*m*cmwagem[m-1]*cmwagem[m+1]
}
cmwagem <- cmwagem[1:4]
secmwage <- sqrt(secmwage/n)
gmean <- exp(mlwage)
segmean <- gmean*slwage/sqrt(n)
ratio <- mwage/meducation
vr <- (swage/meducation)^2 - 2*c[2,1]*mwage/(meducation^3) + ((seducation*mwage)^2)/(meducation^4)
ser <- sqrt(vr/n)
sdwage2 <- sd((wage-mwage)^2)
sevwage <- sdwage2/sqrt(n)
# Gamma Model
alpha <- (mwage^2)/vwage
beta <- vwage/mwage
H <- matrix(0,2,2)
H[1,1] <- -1/beta^2
H[2,1] <- 2*(1+alpha)/beta
H[1,2] <- 1/mwage
H[2,2] <- -(2+1/alpha)
V <- matrix(0,2,2)
V[1,1] <- var(wage^2)
V[2,1] <- V[1,2] <- cov(wage,wage^2)
V[2,2] <- var(wage)
Vg <- t(H)%*%V%*%H/n
# Lognormal
theta <- mlwage
v2 <- vlwage
# student t
mu4 <- mean((wage-mwage)^4)
v4 <- vwage^2
thetat <- mwage
v2t <- mu4*vwage/(2*mu4-3*v4)
r <- (4*mu4-6*v4)/(mu4-3*v4)
length(secmwage)
4
Page 226. Sample mean, covariance matrix and standard errors#
cat("\n")
cat("Sample size, mean wage, mean education \n")
print(c(n,mwage,meducation))
Sample size, mean wage, mean education
[1] 529.00000 31.91416 13.96597
cat("Variance matrix: wage, education \n")
print(c)
cat("\n")
Variance matrix: wage, education
[,1] [,2]
[1,] 834.22452 34.553063
[2,] 34.55306 7.165507
cat("s.e. of mean estimator: wage, education \n")
print(c(sewage,seeducation))
cat("\n")
s.e. of mean estimator: wage, education
[1] 1.2557802 0.1163846
Page 227. Table 11.1#
cat("moments, 1-4, with standard errors \n")
print(mwagem)
print(semwagem)
cat("\n")
moments, 1-4, with standard errors
[1] 0.3191416 0.1851161 0.1868959 0.2539828
[1] 0.01255780 0.02039918 0.03345857 0.05625894
cat("centered moments, 1-4, with standard errors \n")
print(cmwagem)
print(secmwage)
cat("\n")
centered moments, 1-4, with standard errors
[1] 0.31914162 0.08326475 0.07467113 0.09740267
[,1] [,2] [,3] [,4]
[1,] 0.0125578 0.01308984 0.01522651 0.02176415
Page 228#
cat("Geometric mean, sd, & standard error \n")
print(c(gmean,slwage,segmean))
cat("\n")
Geometric mean, sd, & standard error
[1] 24.9376887 0.6627225 0.7185551
Page 229#
cat("ratio of means, variance, & standard error \n")
print(c(ratio,vr,ser))
cat("\n")
ratio of means, variance, & standard error
[1] 2.28513692 3.65921844 0.08316993
cat("variance: variance of X^2, s.e. \n")
print(c(sdwage2,sevwage))
cat("\n")
variance: variance of X^2, s.e.
[1] 3010.6643 130.8984
Page 235#
cat("Gamma Model: alpha, beta \n")
print(c(alpha,beta))
print(sqrt(diag(Vg)))
cat("\n")
Gamma Model: alpha, beta
[1] 1.220911 26.139634
[1] 0.1194149 3.2591115
Page 236#
cat("Lognormal Model: theta, v2 \n")
print(c(theta,v2))
cat("\n")
Lognormal Model: theta, v2
[1] 3.2163803 0.4392011
cat("student t Model: theta, v2, r \n")
print(c(thetat,v2t,r))
cat("\n")
student t Model: theta, v2, r
[1] 31.914162 467.181618 4.545651
Figure 11.1 Best Unbiased Estimation#
suppressPackageStartupMessages(library(spatstat))
# True and Auxiliary Densities
x <- seq(-1,1,.01)
f1 <- (1 - (x^2))*3/4
f2 <- f1*(1+x)
leg1 <- expression(f(x))
leg2 <- expression(f[mu](x))
wd <- 1.4
plot(x,f1,type="l",lty=1,xaxs="i",yaxs="i",xaxt="n",xlim=c(-1,1),ylim=c(-0.075,1),ylab="",xlab="",yaxt="n",bty="n",lwd=wd)
abline(h=0,lwd=wd)
lines(x,f2,lwd=wd)
arrows(0,.75,0,0,angle=20,length=.2)
arrows(.2,108/125,.2,0,angle=20,length=.2)
text(0,-.05,expression(mu[0]))
text(.2,-.05,expression(mu))
text(-.65,.6,leg1)
text(.85,.7,leg2)
# Space of Distribution Functions
x <- seq(-1,1,.01)
y <- (abs(x)^(1/2))*((x>0) - (x<0))
plot.new()
plot.window(c(-4,4),c(-4,4))
W <- ellipse(a=3,b=3.5,centre=c(0,0),phi=0,npoly=1024)
plot(W,add=TRUE,lwd=wd)
lines(x,y,lwd=wd)
points(0,0,pch=19,col="black")
text(-2.2,3.2,expression(bold(F[2])))
text(0.5,0,expression(italic(F[0])))
text(1.4,1,expression(italic(F[mu])))
Figure 11.2. ParametricModels fit to Wage Data using Method of Moments#
dat <- read.table("cps09mar.txt")
experience <- dat[,1]-dat[,4]-6
mwm <- (dat[,11]==1)&(dat[,12]<=2)&(dat[,2]==0)&(experience==20)
dat1 <- dat[mwm,]
wage <- as.matrix(dat1[,5]/(dat1[,6]*dat1[,7]))
x <- as.vector(seq(0.1,100,.1))
# Normal
m <- as.vector(mean(wage))
s <- as.vector(sd(wage))
f1 <- dnorm(x,m,s)
# Gamma
v <- as.vector(var(wage))
alpha <- (m^2)/v
beta <- v/m
f2 <- (x^(alpha-1))*exp(-x/beta)/gamma(alpha)/(beta^alpha)
# Log Normal
lwage <- log(wage)
m2 <- as.vector(mean(lwage))
v2 <- as.vector(var(lwage))
f3 <- exp(-((log(x)-m2)^2)/2/v2)/x/sqrt(2*pi*v2)
# Student t
mu4 <- as.vector(mean((wage-m)^4))
v4 <- v^2
v2t <- mu4*v/(2*mu4-3*v4)
r <- (4*mu4-6*v4)/(mu4-3*v4)
r1 <- (r+1)/2
f4 <- ((((x-m)^2)/r/v2t+1)^(-r1))*gamma(r1)/sqrt(r*pi*v2t)/gamma(r/2)
leg1 <- "Normal"
leg2 <- "Gamma"
leg3 <- "Log Normal"
leg4 <- "Student t"
plot(x,f1,type="l",lty=1,xaxs="i",yaxs="i",ylab="",xlab="",xlim=c(-1,100),ylim=c(0,.031),xaxt="n",yaxt="n",bty="n")
lines(x,f2,lty=5)
lines(x,f3,lty=2)
lines(x,f4,lty=6)
legend("topright",legend=c(leg1,leg2,leg3,leg4),lty=c(1,5,2,6),y.intersp=1.3,cex=.75,bty="n")
axis(side=1,seq(0,100,10),cex.axis=.75)
abline(v=0)
Figure 11.3: Empirical Distribution and Quantile Functions#
dat <- read.table("cps09mar.txt")
experience <- dat[,1]-dat[,4]-6
mbf <- (dat[,11]==2)&(dat[,12]<=2)&(dat[,2]==1)&(experience==12)
dat1 <- dat[mbf,]
wage <- as.matrix(dat1[,5]/(dat1[,6]*dat1[,7]))
w <- sort(wage)
n <- length(wage)
wd <- 1.4
# Empirical Distribution Function
plot.new()
plot.window(c(3,60),c(0.04,1))
points(w[1],0,pch=1,col="black",cex=.6)
points(w[n],1,pch=19,col="black",cex=.6)
for (i in 1:(n-1)) {
points(w[i+1],i/n,pch=1,col="black",cex=.6)
points(w[i],i/n,pch=19,col="black",cex=.6)
lines(c(w[i],w[i+1]-.25),c(i/n,i/n),lwd=wd)
points(w[i],0.02,pch=4,col="black",cex=.6)
}
points(w[n],0.02,pch=4,col="black",cex=.6)
lines(c(0,w[1]-.2),c(0,0),lwd=wd)
lines(c(w[n],60),c(1,1),lwd=wd)
axis(side=1,seq(-10,70,10),lwd=wd)
axis(side=2,seq(-1,1,.2),lwd=wd)
# Empirical Quantile Function
plot.new()
plot.window(c(.03,.97),c(3,60))
abline(h=0,lwd=wd)
for (i in 1:(n)) {
points((i-1)/n,w[i],pch=19,col="black",cex=.6)
points((i)/n,w[i],pch=1,col="black",cex=.6)
lines(c((i-1)/n,(i)/n-.005),c(w[i],w[i]),lwd=wd)
}
points(1,w[n],pch=19,col="black",cex=.6)
axis(side=1,seq(-1,1,.2),lwd=wd)
axis(side=2,seq(-10,60,10),lwd=wd)
