Quick visual analysis construction vs demolition activity, per neighborhood:
my.plot.theme <- theme(text = element_text(family = "Futura Bk BT"),
#plot.title = element_text(size = 36),
axis.line = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
#axis.text.y = element_blank(),
#legend.title = element_text(size = 16),
#legend.text = element_text(size = 24),
#strip.text = element_text(size = 20),
panel.grid.major = element_blank (), # remove major grid
panel.grid.minor = element_blank (), # remove minor grid
axis.text = element_blank (),
axis.title = element_blank (),
axis.ticks = element_blank ())
ggplot(data = BPdf_by_NB) +
geom_line(aes(x = YEARMONTH, y = NEWCONS, group = 1, color = "new construction")) +
geom_line(aes(x = YEARMONTH, y = DEMOLITIONS, group = 1, color = "demolition")) +
scale_colour_manual(values=c("new construction"="darkseagreen3","demolition" = "black"),
guide = guide_legend(title = NULL)) +
theme_minimal()+
theme (legend.position = "bottom") +
labs(title="New constructions Vs Demolitions (2010 - 2015)\n") +
my.plot.theme +
facet_wrap(~ NB)
We will concentrate on demolitions in the Downtown (Central district), as the comparative plot by neighboorhood shows a much higher occurrence there than anywhere else in the city.
Downtown.demolitions.ts <- ts(BPdf_by_NB[BPdf_by_NB$NB == "Central", "DEMOLITIONS"], start = c(2010,1), frequency = 12)
hchart(Downtown.demolitions.ts) %>%
hc_title(text = "Approved Demolition permits - Boston Downtown")
Let’s take a look at its first-order differencing:
hchart(diff(Downtown.demolitions.ts))
That took care of the trend. The mean is clearlt stable now, but there’s a worrying increase of the variance in later times, that may indicate non-stationarity. To be sure, we will take a root’unit test to check for non-stationarity:
tseries::adf.test(diff(Downtown.demolitions.ts), k = 12)
##
## Augmented Dickey-Fuller Test
##
## data: diff(Downtown.demolitions.ts)
## Dickey-Fuller = -3.9278, Lag order = 12, p-value = 0.0184
## alternative hypothesis: stationary
With a p-value well below the threshold at 0.5, we can assume that the differenced time series is stationary. We can continue with a look at the autocorrelation plots:
tsdisplay(diff(Downtown.demolitions.ts))
There’s significant autocorrelation at lag 1, suggesting a MA(1) process. There is also significant autocorrelation at lags 14 and 15, perhaps as a result of seasonality. There’s partial autocorrelation at lags 1 and 2 -AR(2)- and spike of partial autocorrelation at lag 13, indicating the possibility of yearly AR seasonality.
AR (autoregression) models predicts future values of a variable of interest using a linear combination of past values of the variable. The term “autoregression” is used to indicate that it is a regression of the variable against itself.
An autoregressive model of order \(p\) is written as:
\[ y_{t} = c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + e_{t}, \]
where \(c\) is a constant and \(e_t\) is white noise. This is referred to as an AR(\(p\)) model.
MA (moving average) models, instead of using past values of the forecast variable in a regression, use past forecast errors in a regression-like model.
\[ y_{t} = c + e_t + \theta_{1}e_{t-1} + \theta_{2}e_{t-2} + \dots + \theta_{q}e_{t-q}, \]
where \(e_t\) is white noise. We refer to this as an MA(\(q\)) model.
When a process is modelled with an AR component, and also an MA one, it’s called and ARMA model. When integration (time series differencing) is added, it’s called an ARIMA model.
ARIMA models can also model seasonal data. A seasonal ARIMA model is formed by including additional seasonal terms like this:
Where \(m\) = the number of periods per season. We use uppercase notation for the seasonal parts of the model, and lowercase notation for the non-seasonal parts of the model.
The seasonal part of the model consists of terms that are very similar to the non-seasonal components of the model, but they involve backshifts of the seasonal period. For example, an ARIMA(1,1,1)(1,1,1) model (without a constant) is for quarterly data (\(m=4\)) and can be written as 
(See https://www.otexts.org/fpp/8/9 )
We’ll model the process as an ARIMA(2,1,1)(1,1,0). That is, and AR(2), I(1), MA(1) process with AR and I seasonality:
dd.fit <- Arima(Downtown.demolitions.ts, order = c(2,1,1), seasonal = c(1,1,0))
Now we’ll examine the distribution of the residuals with a custom function (courtesy of http://www.quantlego.com):
PlotForecastErrors <- function(forecasterrors)
{
mybinsize <- IQR(forecasterrors)/4
mysd <- sd(forecasterrors)
mymin <- min(forecasterrors)-mysd*5
mymax <- max(forecasterrors)+mysd*5
mynorm <- rnorm(10000,mean=0,sd=mysd)
mymin2 <- min(mynorm)
mymax2 <- max(mynorm)
if (mymin2<mymin)
{
mymin<-mymin2
}
if (mymax2>mymax)
{
mymax<-mymax2
}
mybins <- seq(mymin,mymax,mybinsize)
hist(forecasterrors, col="red", freq=FALSE, breaks=mybins)
myhist <- hist(mynorm, plot=FALSE, breaks=mybins)
points(myhist$mids, myhist$density, type="l", col="blue", lwd=2)
}
PlotForecastErrors(dd.fit$residuals)
The histogram shows that the distribution of forecast errors is mostly centered on zero, and normally distributed, it is perfectly plausible that the forecast errors are normally distributed with mean zero and constant variance.
tsdisplay(dd.fit$residuals)
The residuals don’t dislay any autocorrelation, so we can assume that our model provides a good approximation.
This is our two-year forecast for approved demolitions at the City’s downtown:
hchart(forecast(dd.fit, h=24)) %>%
hc_title(text = "Two year forecast: Approved Demolition permits - Boston Downtown")
