data-science-platform/R-code-for-RDataMining-book/cs-HPI-forecast.R

   1 ### R code from vignette source 'cs-HPI-forecast.rnw'
   2
   3 ###################################################
   4 ### code chunk number 1: cs-HPI-forecast.rnw:7-10
   5 ###################################################
   6 # free memory
   7 rm(list = ls())
   8 gc()
   9
  10
  11 ###################################################
  12 ### code chunk number 2: cs-HPI-forecast.rnw:31-44
  13 ###################################################
  14 # import data
  15 filepath <- "./data/"
  16 filename <- "House-index-canberra.csv"
  17 houseIndex <- read.csv(paste(filepath, filename, sep=""), header=FALSE)
  18 names(houseIndex) <- c("date", "index")
  19 n <- nrow(houseIndex)
  20 # check start date and end date
  21 cat(paste("HPI from", houseIndex$date[1], "to", houseIndex$date[n], "\n"))
  22 # extract year and month
  23 dates <- strptime(houseIndex$date, format="%d-%b-%y")
  24 houseIndex$year <- dates$year + 1900
  25 houseIndex$month <- dates$mon + 1
  26 fromYear <- houseIndex$year[1]
  27
  28
  29 ###################################################
  30 ### code chunk number 3: cs-HPI-forecast.rnw:48-51 (eval = FALSE)
  31 ###################################################
  32 ## dates <- as.Date(houseIndex$date, format="%d-%b-%y")
  33 ## houseIndex$year <- as.numeric(format(dates, "%Y"))
  34 ## houseIndex$month <- as.numeric(format(dates, "%m"))
  35
  36
  37 ###################################################
  38 ### code chunk number 4: cs-HPI-forecast.rnw:63-74
  39 ###################################################
  40 plot(houseIndex$index, pty=1, type="l", lty="solid", xaxt="n", xlab="",
  41      ylab="Index", main=paste("HPI (Canberra) - Since ", fromYear, sep=""))
  42 # draw tick-marks at 31 Jan of every year
  43 nYear <- ceiling(n/12)
  44 posEveryYear <- 12 * (1:nYear) - 11
  45 axis(1, labels=houseIndex$date[posEveryYear], las=3, at=posEveryYear)
  46 # add horizontal reference lines
  47 abline(h=1:4, col="gray", lty="dotted")
  48 # draw a vertical reference line every five years
  49 posEvery5years <- 12 * (5* 1:ceiling(nYear/5) - 4) - 11
  50 abline(v=posEvery5years, col="gray", lty="dotted")
  51
  52
  53 ###################################################
  54 ### code chunk number 5: cs-HPI-forecast.rnw:88-93
  55 ###################################################
  56 houseIndex$delta <- houseIndex$index - c(1, houseIndex$index[-n])
  57 plot(houseIndex$delta, main="Increase in HPI", xaxt="n", xlab="")
  58 axis(1, labels=houseIndex$date[posEveryYear], las=3, at=posEveryYear)
  59 # add a reference line
  60 abline(h=0, lty="dotted")
  61
  62
  63 ###################################################
  64 ### code chunk number 6: cs-HPI-forecast.rnw:109-119
  65 ###################################################
  66 # increase ratio in every month
  67 houseIndex$rate <- houseIndex$index / c(1, houseIndex$index[-n]) - 1
  68 # percentage of months having positive increases in HPI
  69 100 * sum(houseIndex$rate>0)/n
  70 # use ifelse() to set positive values to green and and negative ones to red
  71 plot(houseIndex$rate, xaxt="n", xlab="", ylab="HPI Increase Rate",
  72      col=ifelse(houseIndex$rate>0,"green","red"),
  73      pch=ifelse(houseIndex$rate>0,"+","o"))
  74 axis(1, labels=houseIndex$date[posEveryYear], las=3, at=posEveryYear)
  75 abline(h=0, lty="dotted")
  76
  77
  78 ###################################################
  79 ### code chunk number 7: cs-HPI-forecast.rnw:134-141
  80 ###################################################
  81 rateMatrix <- xtabs(rate ~ month + year, data=houseIndex)
  82 # show the first four years, rounded to 4 decimal places
  83 round(rateMatrix[,1:4], digits=4)
  84 # plot a grouped barchart:
  85 barplot(rateMatrix, beside=TRUE, space=c(0,2),
  86         col=ifelse(rateMatrix>0,"lightgreen","lightpink"),
  87         ylab="HPI Increase Rate", cex.names=1.2)
  88
  89
  90 ###################################################
  91 ### code chunk number 8: cs-HPI-forecast.rnw:153-155
  92 ###################################################
  93 numPositiveMonths <- colSums(rateMatrix>0)
  94 barplot(numPositiveMonths, xlab="Year", ylab="Number of Months with Increased HPI")
  95
  96
  97 ###################################################
  98 ### code chunk number 9: cs-HPI-forecast.rnw:164-167
  99 ###################################################
 100 yearlyMean <- colMeans(rateMatrix)
 101 barplot(yearlyMean, main="Yearly Average Increase Rates of HPI",
 102         col=ifelse(yearlyMean>0,"lightgreen","lightpink"), xlab="Year")
 103
 104
 105 ###################################################
 106 ### code chunk number 10: cs-HPI-forecast.rnw:178-181
 107 ###################################################
 108 monthlyMean <- rowMeans(rateMatrix)
 109 plot(names(monthlyMean), monthlyMean, type="b", xlab="Month",
 110      main="Monthly Average Increase Rates of HPI")
 111
 112
 113 ###################################################
 114 ### code chunk number 11: cs-HPI-forecast.rnw:193-195
 115 ###################################################
 116 summary(houseIndex$rate)
 117 boxplot(houseIndex$rate, ylab="HPI Increase Rate")
 118
 119
 120 ###################################################
 121 ### code chunk number 12: cs-HPI-forecast.rnw:205-206
 122 ###################################################
 123 boxplot(rate ~ year, data=houseIndex, xlab="Year", ylab="HPI Increase Rate")
 124
 125
 126 ###################################################
 127 ### code chunk number 13: cs-HPI-forecast.rnw:215-216
 128 ###################################################
 129 boxplot(rate ~ month, data=houseIndex, xlab="Month", ylab="HPI Increase Rate")
 130
 131
 132 ###################################################
 133 ### code chunk number 14: cs-HPI-forecast.rnw:234-237
 134 ###################################################
 135 hpi <- ts(houseIndex$index, start=c(1990,1), frequency=12)
 136 f <- stl(hpi, "per")
 137 plot(f)
 138
 139
 140 ###################################################
 141 ### code chunk number 15: cs-HPI-forecast.rnw:247-250
 142 ###################################################
 143 # plot seasonal components
 144 plot(f$time.series[1:12,"seasonal"], type='b', xlab="Month",
 145      ylab="Seasonal Components")
 146
 147
 148 ###################################################
 149 ### code chunk number 16: cs-HPI-forecast.rnw:258-263 (eval = FALSE)
 150 ###################################################
 151 ## # an alternative decomposition function
 152 ## f2 <- decompose(hpi)
 153 ## plot(f2)
 154 ## # plot seasonal components
 155 ## plot(f2$figure, type="b", xlab="Month", ylab="Seasonal Components")
 156
 157
 158 ###################################################
 159 ### code chunk number 17: cs-HPI-forecast.rnw:277-298
 160 ###################################################
 161 startYear <- 1990
 162 endYear <- 2010
 163 # to forecast HPIs in the next four years
 164 nYearAhead <- 4
 165
 166 fit <- arima(hpi, order=c(2,0,1), seasonal=list(order=c(2,1,0), period=12))
 167 fore <- predict(fit, n.ahead=12*nYearAhead)
 168 # error bounds at 95% confidence level
 169 U <- fore$pred + 2 * fore$se
 170 L <- fore$pred - 2 * fore$se
 171 # plot original and predicted data, as well as error bounds
 172 ts.plot(hpi, fore$pred, U, L, col=c("black", "blue","green","red"),
 173         lty=c(1,5,2,2), gpars=list(xaxt="n",xlab=""),
 174         ylab="Index", main="House Price Trading Index Forecast (Canberra)")
 175 # add labels, reference grid and legend
 176 years <- startYear:(endYear+nYearAhead+1)
 177 axis(1, labels=paste("Jan ", years, sep=""), las=3, at=years)
 178 grid()
 179 legend("topleft", col=c("black", "blue","green","red"), lty=c(1,5,2,2),
 180        c("Actual Index", "Forecast", "Upper Bound (95% Confidence)",
 181          "Lower Bound (95% Confidence)"))
 182
 183
 184 ###################################################
 185 ### code chunk number 18: cs-HPI-forecast.rnw:308-317
 186 ###################################################
 187 ts.plot(fore$pred, U, L, col=c("blue","green","red"),
 188         lty=c(5,2,2), gpars=list(xaxt="n",xlab=""),
 189         ylab="Index", main="House Price Trading Index Forecast (Canberra)")
 190 years <- endYear + (1 : (nYearAhead+1))
 191 axis(1, labels=paste("Jan ", years, sep=""), las=3, at=years)
 192 grid(col = "gray", lty = "dotted")
 193 legend("topleft", col=c("blue","green","red"), lty=c(5,2,2),
 194        c("Forecast", "Upper Bound (95% Confidence)",
 195          "Lower Bound (95% Confidence)"))
 196
 197
 198 ###################################################
 199 ### code chunk number 19: cs-HPI-forecast.rnw:327-335
 200 ###################################################
 201 newHpi <- ts(c(hpi, fore$pred), start=c(1990,1), frequency=12)
 202 (startDate <- start(newHpi))
 203 startYear <- startDate[1]
 204 m <- 9 + (2009-startYear)*12
 205 n <- 9 + (2011-startYear)*12
 206 # percentage of increase
 207 100 * (newHpi[n] / newHpi[m] - 1)
 208 round(535000 * newHpi[n] / newHpi[m])
 209
 210