# Example 10.9, Puerto Rican Employment
# Data set: prminwge

# Function for result reporting
source("_report.R")

# Load the data and estimate the model in the background
load("prminwge.Rdata")
names(data)[names(data)=="lprepop"]="lprepop.t"
names(data)[names(data)=="lmincov"]="lmincov.t"
names(data)[names(data)=="lusgnp"]="lusgnp.t"
data=ts(data,start=1950,frequency=1)
model0=lm(lprepop.t~lmincov.t+lusgnp.t,data=data)
model1=lm(lprepop.t~lmincov.t+lusgnp.t+t,data=data)
dig0=c(2,3,3,3)
dig1=c(2,3,2,3,3)

# Recap model
{
  cat("This example closely follows Example 10.3. The estimation results in the previous example were")
  reportreg(model0,dig0,suffix=".hat",adj=T)
}

# Estimate first model
{
  cat("Adding a linear time trend into the model gives")
  reportreg(model1,dig1,suffix=".hat",adj=T)
}

# Interpretation
m=lm(lusgnp.t~t,data=data)
cat("The coefficient on lusgnp has changed from ", printcoef(model0,3,dig0[3]), " and insignificant to ",
    printcoef(model1,3,dig1[3]), " and very significant. The change in the coefficient on lmincov is less notable, but it does become more significant",
    "\nWhile lprepop has no clear upward or downward trend, lusgnp does have a upward linear trend: regressing lusgnp on t gives a slope coefficient of ",
    printcoef(m,2,2), ", so usgnp grew by about ", 100*as.numeric(printcoef(m,2,2)), "% per year over the period. The coefficient on lusgnp has a detrended interpretation: when usgnp increases by 1% ABOVE its long-run trend, prepop is predicted to increase by ",
    printcoef(model1,3,dig1[3]), "%",
    sep="")