# HINT: # The program assumes that all files are stored in # C:\monographregression\computercode # change directory if code is located elsewhere # to store the figures a subfolder `figures' has to be # created setwd("c:/monographregression/computercode") library("foreign") library("quantreg") library("gamlss") # read data rents <- read.dta("data/stata/rent99.dta") # estimate model for location and scale with linear effects gamlss.rent <- gamlss(formula=rent~area+yearc, sigma.formula=~area+yearc, data=rents, family=NO()) # summarize results summary(gamlss.rent) # estimate model for location and scale with additive effects # preliminary model for smoothing parameter choice # see Lancaster-Booklet.pdf gamlss.rent.add <- quote(gamlss(formula=rent~pb(area,inter=10,degree=3,order=2,df=p[1])+ pb(yearc,inter=10,degree=3,order=2,df=p[2]), sigma.formula=~pb(area,inter=10,degree=3,order=2,df=p[3])+ pb(yearc,inter=10,degree=3,order=2,df=p[4]) , data=rents, family=NO(), trace=FALSE)) # obtain optimal smoothing parameters op.rent <- find.hyper(model=gamlss.rent.add, par=c(3,3,3,3),lower=c(1,1,1,1),upper=c(6,6,6,6), steps=c(0.2,0.2,0.2,0.2),pen=2,trace=FALSE) # final additive model gamlss.rent.add.final <- gamlss(formula=rent~pb(area,inter=10,degree=3,order=2,df=op.rent$par[1])+ pb(yearc,inter=10,degree=3,order=2,df=op.rent$par[2]), sigma.formula=~pb(area,inter=10,degree=3,order=2,df=op.rent$par[3])+ pb(yearc,inter=10,degree=3,order=2,df=op.rent$par[4]) , data=rents, family=NO(), trace=FALSE) # plot estimated functions term.plot(gamlss.rent.add.final,what=c("mu"),se=TRUE) term.plot(gamlss.rent.add.final,what=c("sigma"),se=TRUE) # plot estimated functions and write as postscript files postscript("figures/gamlss_area_mu.eps", height=4.5, width=4.5, horizontal=F) term.plot(gamlss.rent.add.final,what=c("mu"),se=TRUE,terms=1,ylab=" ", xlab="area in sqm",col.se=1,col.term=1,main="mean: area") dev.off() postscript("figures/gamlss_yearc_mu.eps", height=4.5, width=4.5, horizontal=F) term.plot(gamlss.rent.add.final,what=c("mu"),se=TRUE,terms=2, ylab=" ",xlab="year of construction",col.se=1,col.term=1,main="mean: year of construction") dev.off() postscript("figures/gamlss_area_sigma.eps", height=4.5, width=4.5, horizontal=F) term.plot(gamlss.rent.add.final,what=c("sigma"),se=TRUE,terms=1, ylab=" ",xlab="area in sqm",col.se=1,col.term=1,main="standard dev.: area") dev.off() postscript("figures/gamlss_yearc_sigma.eps", height=4.5, width=4.5, horizontal=F) term.plot(gamlss.rent.add.final,what=c("sigma"),se=TRUE,terms=2, ylab=" ",xlab="year of construction",col.se=1,col.term=1,main="standard dev.: year of construction") dev.off()