# R code for the figures for lecture 6

require(np)
require(lattice)
data(oecdpanel)

############################# Figure 1 ############################

# Fit the joint density of logarithmic population growth and investment
# rates
popinv <- npudens(~popgro+inv, data=oecdpanel)
# Call the np plotting routine, but just ask it to calculate, rather than to
# make a plot --- see help(npplot)
fhat <- plot(popinv,plot.behavior="data")
# Extract the one joint density (since it could return multiple objects
# together
fhat <- fhat$d1
  # Now fhat$eval contains the grid of evaluation points, and fhat$dens the
  # actual density values
# Make a contour plot of the joint density, with more levels, and smaller
# labels on the levels, than usual
contourplot(fhat$dens~fhat$eval$Var1*fhat$eval$Var2,cuts=20,xlab="popgro",
            ylab="inv",labels=list(cex=0.5))


############################# Figure 2 ############################

# Fit the density of logarithmic population growth rates conditional on year
pop.cdens <- npcdens(popgro ~ year,data=oecdpanel)
# Manually construct a grid of points on which to plot
plotting.grid <- expand.grid(year=seq(from=1965,to=1995,by=1),
                             popgro=seq(from=-3.5,to=-2.4,length.out=300))
# Evaluate the function
fhat <- predict(pop.cdens,newdata=plotting.grid)
# Make a wireframe plot
wirefame(fhat~plotting.grid$year*plotting.grid$popgro,scales=list(arrows=FALSE),
         xlab="year",ylab="popgro",zlab="pdf")



############################ Figure 3 ##############################

# Fit the density of popgro conditional on year and OECD membership
pop.cdens.o <- npcdens(popgro~year+factor(oecd),data=oecdpanel)
  # Ensure that npcdens treats oecd as a categorical variable, not a number
# Make the grid, but now include the OECD factor
plotting.grid <- expand.grid(year=seq(from=1965,to=1995,by=1),
                               popgro=seq(from=-3.4,to=-2.4,length.out=300),
                               oecd=unique(oecdpanel$oecd))
fhat <- predict(pop.cdens.o,newdata=plotting.grid)
# Make side-by-side wireframes, conditional on OECD membership
wireframe(fhat~plotting.grid$year*plotting.grid$popgro|plotting.grid$oecd,
          scales=list(arrows=FALSE),xlab="year",ylab="popgro",zlab="pdf")


############################ Figure 4 ##############################

# Code omitted --- reproducing the figure is a challenge exercise