if (FALSE) {
## This illustration was made possible by Samuele Centorrino
##
set.seed(42)
n <- 1500
## The DGP is as follows:
## 1) y = phi(z) + u
## 2) E(u|z) != 0 (endogeneity present)
## 3) Suppose there exists an instrument w such that z = f(w) + v and
## E(u|w) = 0
## 4) We generate v, w, and generate u such that u and z are
## correlated. To achieve this we express u as a function of v (i.e. u =
## gamma v + eps)
v <- rnorm(n,mean=0,sd=0.27)
eps <- rnorm(n,mean=0,sd=0.05)
u <- -0.5*v + eps
w <- rnorm(n,mean=0,sd=1)
## In Darolles et al (2011) there exist two DGPs. The first is
## phi(z)=z^2 and the second is phi(z)=exp(-abs(z)) (which is
## discontinuous and has a kink at zero).
fun1 <- function(z) { z^2 }
fun2 <- function(z) { exp(-abs(z)) }
z <- 0.2*w + v
## Generate two y vectors for each function.
y1 <- fun1(z) + u
y2 <- fun2(z) + u
## You set y to be either y1 or y2 (ditto for phi) depending on which
## DGP you are considering:
y <- y1
phi <- fun1
## Create an evaluation dataset sorting on z (for plotting)
evaldata <- data.frame(y,z,w)
evaldata <- evaldata[order(evaldata$z),]
## Compute the non-IV regression spline estimator of E(y|z)
model.noniv <- crs(y~z,opts=opts)
mean.noniv <- predict(model.noniv,newdata=evaldata)
## Compute the IV-regression spline estimator of phi(z)
model.iv <- crsiv(y=y,z=z,w=w)
phi.iv <- predict(model.iv,newdata=evaldata)
## For the plots, restrict focal attention to the bulk of the data
## (i.e. for the plotting area trim out 1/4 of one percent from each
## tail of y and z)
trim <- 0.0025
curve(phi,min(z),max(z),
xlim=quantile(z,c(trim,1-trim)),
ylim=quantile(y,c(trim,1-trim)),
ylab="Y",
xlab="Z",
main="Nonparametric Instrumental Spline Regression",
sub=paste("Landweber-Fridman: iterations = ", model.iv$num.iterations,sep=""),
lwd=1,lty=1)
points(z,y,type="p",cex=.25,col="grey")
lines(evaldata$z,evaldata$z^2 -0.325*evaldata$z,lwd=1,lty=1)
lines(evaldata$z,phi.iv,col="blue",lwd=2,lty=2)
lines(evaldata$z,mean.noniv,col="red",lwd=2,lty=4)
legend(quantile(z,trim),quantile(y,1-trim),
c(expression(paste(varphi(z),", E(y|z)",sep="")),
expression(paste("Nonparametric ",hat(varphi)(z))),
"Nonparametric E(y|z)"),
lty=c(1,2,4),
col=c("black","blue","red"),
lwd=c(1,2,2))
}
Run the code above in your browser using DataLab