An Evaluation of 4 Ways to Forecast Daily Time Series with Annual Seasonality
Clayton Glasser
March 15, 2019
The purpose of this experiment is to compare 4 methods of creating a 1-year-out forecast of precipitation levels based on data that consists of daily measurements of precipitation volume over a 10-year period from a single location. Additional dimensions of data are assumed unavailable, so multiple linear regression based techniques are not an option. However each measurement is associated with a date; therefore, it can can be treated as a univariate time series.
Time series are often found in econometrics and finance, where their dynamics can be very different than those found in a hydrological context, where calendar dates and cultural phenomena have little influence. Numerous methods exist for analysing and forecasting time series. The ones chosen for this experiment were chosen partially on the basis on their ability to handle the unique characteristics of this precipitation data – particularly the long seasonal period of 365 days, which exceeds the computational capacity of some common time series related functions, which were designed with the smaller frequencies relevant to cultural data in mind, e.g. monthly (12) and quarterly (4).
The precipitation data contains 365 measurements per cycle repetition (technically, the average length of a year is 365.25, because of leap years, but it is treated as 356 here). The primary seasonal cycle associated with this precipitation data is an annual cycle correlating to the intensification of rainfall which zeniths in winter and nadirs in summer.
The four methods of analysis and forecasting evaluated here are:
- Seasonal ARIMA
- Holt-Winters Seasonal Exponential Smoothing
- Non-seasonal ARIMA with a Fourier transform as an exogenous covariate
- Christiano-Fitzgerald finite data approximation to the ideal time series bandpass filter
The results of these techniques are evaluated relative to each other by measuring the root mean squared error (RMSE) of 1-year forecasts of the 10th year of a data (out-of-sample) based on training sets of data consisting of the first 9 years of data, as well as of the training data itself (in-sample).
Data Preparation
Required Packages
This project requires the following R packages: multiple packages found in the tidyverse package, including dplyr (piping), ggplot2 (plotting), and readr (file reading); forecast is used for time series creation and forecasting; mfilter for the creation of band pass filters; lubridate for date management; tseries for testing functions.
library(tidyverse)
library(forecast)
library(mFilter)
library(lubridate)
library(tseries)Read in the Data
The data used here has already been tidied up and wrangled into an efficient structure. X1 is an index. Date contains the calendar date of the measurement. water_date contains the equivalent hydrological date to the calendar date (rainfall chronology runs Oct-Sep). PRCP contains the precipitation volume measurement.
rainfall <- read_csv("US1ORLA0076.csv")## Parsed with column specification:
## cols(
## X1 = col_integer(),
## DATE = col_date(format = ""),
## water_date = col_date(format = ""),
## PRCP = col_double()
## )head(rainfall)## # A tibble: 6 x 4
## X1 DATE water_date PRCP
## <int> <date> <date> <dbl>
## 1 1 2008-10-01 2007-10-01 0
## 2 2 2008-10-02 2007-10-02 0.02
## 3 3 2008-10-03 2007-10-03 0.31
## 4 4 2008-10-04 2007-10-04 1.8
## 5 5 2008-10-05 2007-10-05 0.03
## 6 6 2008-10-06 2007-10-06 0.19Visualize the raw data
ggplot(rainfall, aes(water_date, PRCP)) +
geom_point() +
scale_x_date(date_labels = "%b %y", date_breaks = "6 months", date_minor_breaks = "3 months") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
A few interesting qualities are clear upon visual inspection.
Because rainfall frequently does not occur, but is always additive when it does, many of the data points are 0 or very close to 0, but never go below 0.
There is a clear presence of at least one seasonal cycle corresponding to the length of one year, repeated 10 times, once for each of the 10 years. Intra-annual cycles may exist but are not visibly apparent.
The data appears to exhibit low-frequency characteristics, visible in its modulating sharply between large and small values.
Data Transformation
The following techniques are predicated on the understanding of our data as comprised of “signals” and “noise”. The signal components are abstractions of the literal data into simplified shapes that represent it’s essential form or constituent forms. In this case, the signals will be linear regressions of the data’s long-term trend and short-term cycles. The “remainder” component is every aspect of the data that is not represented by these linear regressions. If the signals are perfectly captured, the remainder will be perfectly Gaussian noise; if they are not, some signal will remain in the remainder and it will not be perfectly Gaussian.
To aid in the process of identifying the constituent signals, we first apply two forms of non-destructive and reversible transformation to the data, which compresses it down a less complex expression, with a more essential shape, which will improve model selection.
Standardization
Subtract the mean and divide by the standard deviation.
rainfall$prcp_scaled <- scale(rainfall$PRCP)Normalization
We use the inverse hyperbolic sin function as the preferred method to log-transform data with many zeros, which this data has because rainfall has a natural valence of zero.
rainfall$prcp_log <- log(rainfall$prcp_scaled + sqrt(rainfall$prcp_scaled ^ 2 + 1))Visualize the transformed data
After standardization, the data is essentially z-scores, and the scale no longer represents the absolute quantity of precipitation (blue), but the value of the data’s deviation from the mean (black). Thus, some measurements register below zero even though sub-zero quantities of rainfall are not possible.
ggplot() +
geom_point(data=rainfall, aes(water_date, prcp_log, colour="Transformed"),alpha = 0.33) +
geom_point(data=rainfall, aes(water_date, PRCP, colour="Original"), alpha = 0.33) +
scale_x_date(date_labels = "%b %y", date_breaks = "6 months", date_minor_breaks = "3 months") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_colour_manual(name="Data Iteration", values=c(Original="blue", Transformed="black")) +
ylab("Precipitation") + xlab("Time")
Create time series object
Most of the analyses and forecasting methods used below require time series objects as inputs, so we construct one using the transformed data which spans a precisely 10 year long period, with a frequency of 365, as a close approximation of the technically accurate 365.25 period.
rainfall_transformed_ts <- ts(rainfall$prcp_log, start=c(2007, 10), end=c(2017, 9), frequency = 365)
plot(rainfall_transformed_ts)
Forecasting with ARIMA
ARIMA models are the most general class of models for forecasting those time series which can be made to be “stationary” by differencing. A stationary series is defined as having no trend (long term statistically significant incline or decline) and varying at a constant amplitude around its mean, which means that its autocorrelations (correlations with its own prior deviations from the mean) remain constant over time. Time series of this form can be viewed as a combination of signal and noise, with the signal(s) being a pattern of fast or slow mean reversion (sinusoidal oscillation), or rapid alternation in sign. An ARIMA model can be viewed as a “filter” that tries to separate the signal from the noise, and the signal is then extrapolated into the future to obtain forecasts.
The acronym ARIMA stands for Auto-Regressive Integrated Moving Average. Lags of the stationarized series in the forecasting equation are called “autoregressive” terms (AR), lags of the forecast errors are called “moving average” terms (MA), and a time series which needs to be differenced to be made stationary is said to be an “integrated” version of a stationary series (I). Exponential smoothing models, which are explored in the following section, are special cases of ARIMA models.
The ARIMA forecasting equation for a stationary time series is a linear (i.e., regression-type) equation in which the predictors consist of lags of the dependent variable and/or lags of the forecast errors.
Build ARIMA model
Decompose the Time Series
Decompose the time series into seasonal, trend and irregular (remainder, noise) components using loess, available in the stl() function of forecast.
rainfall_decomp1 = stl(rainfall_transformed_ts, s.window = "periodic")
plot(rainfall_decomp1)
We can see by inspecting the plot of decomposed elements that a clear annual seasonal pattern has been detected, but it seems that some degree of seasonality is visible in the residuals, meaning that the LOESS process has not completely extracted the signal, and has left some of it behind in the remainder, which would ideally be white noise.
It is also interesting to note the substantial dip in the trend component around 2013, which makes sense, as this precipitation data is from coastal Oregon, and 2013 was the driest year on record for large parts of the west coast of the US: https://www.wunderground.com/blog/weatherhistorian/driest-year-on-record-for-california-oregon-wettest-in-asheville-ma.html
Deseasonalize the series by extracting the identified trend and seasonality
rainfall_deseasonal1 <- seasadj(rainfall_decomp1)
plot(rainfall_deseasonal1)
The resulting data is what will be fed into the ARIMA filter. As expected, there appears to be some seasonal information still present in this data, suggesting our prediction will be naive.
Testing for stationarity
The augmented Dickey–Fuller test (ADF) tests the null hypothesis that a unit root is present in a time series sample. The alternative hypothesis is trend-stationarity.
adf.test(rainfall_deseasonal1, alternative = "stationary") ##
## Augmented Dickey-Fuller Test
##
## data: rainfall_deseasonal1
## Dickey-Fuller = -12.917, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationaryThe test concludes that the data is stationary, meaning no differencing will be required for this ARIMA model, or in other words, the I term of the ARIMA model will be of order 0.
Testing for autocorrelation
It is usually not possible to tell from a time plot what values of p (AR) and q (MA) are appropriate for the data. However, it is sometimes possible to use the ACF plot, and the closely related PACF plot, to estimate appropriate values for p and q.
The ACF plot is a bar chart of the coefficients of correlation between a time series and lags of itself. The PACF plot is a plot of the partial correlation coefficients between the series and lags of itself.
Acf(rainfall_deseasonal1, main='',lag.max=1000)
Pacf(rainfall_deseasonal1, main='',lag.max=1000)
Autocorrelation still present, with exponentially diminishing lags appearing at the beginning of the series, suggesting the forecasting may be aided by at least one AR term, and possibly more.
Using auto.arima to automatically select optimal ARIMA p and q terms.
The auto.arima function from the R’s forecast package estimates the ARIMA model using maximum likelihood estimation (MLE) which finds the values of the parameters which maximise the probability of obtaining the data that we have observed. In practice, R will report the value of the log likelihood of the data; that is, the logarithm of the probability of the observed data coming from the estimated model.
rainfall_arima1 <- auto.arima(rainfall_deseasonal1, seasonal=FALSE, approximation = FALSE)
rainfall_arima1 ## Series: rainfall_deseasonal1
## ARIMA(1,0,2) with non-zero mean
##
## Coefficients:
## ar1 ma1 ma2 mean
## 0.7610 -0.4411 -0.0788 -0.0863
## s.e. 0.0397 0.0437 0.0244 0.0178
##
## sigma^2 estimated as 0.2878: log likelihood=-2904.22
## AIC=5818.45 AICc=5818.46 BIC=5849.46auto.arima() asserts an ARIMA model with the expected one AR term, no differencing, and two MA terms: ARIMA(1,0,2). This could also be referred to as an ARMA model, since there is no I (d) term.
If we inspect the residuals present in the ARIMA model, we can get a sense for how much of the signal it has captured.
checkresiduals(rainfall_arima1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2) with non-zero mean
## Q* = 853.69, df = 726, p-value = 0.0007118
##
## Model df: 4. Total lags used: 730Pattern visible in the residuals time plot, lag spikes at the beginning of the frequency in the ACF plot, and the non-Gaussian shape of the parametrically measured residuals all suggest some correlation, or signal, is not captured by the model.
Forecasting using decomposition by LOESS and auto.arima modeling
Plot WITHOUT extracted seasonality
plot(forecast(rainfall_arima1, h=365)) 
Plot WITH extracted seasonality
rainfall_arima1_w_seas <- auto.arima(rainfall_deseasonal1, seasonal=TRUE, approximation = FALSE) #create ARIMA model
plot(forecast(rainfall_arima1_w_seas, h=365), xlab="Date", ylab="Precipitation (z-score)") #plot forecast
The forecast both with and without seasonality incorporated is very naive, apparently just a projection of the mean. I believe this is because the auto.arima() function will allow a seasonal period up to 350 (not 365) but in practice will usually run out of memory when the seasonal period is more than about 200.
Testing forecast made with seasonal decomposition by LOESS and auto.arima
This forecast isn’t going to win any competitions, but we’ll measure the root mean squared error so we can compare it to the following methods.
# Training set - 9 years
rainfall_train1 <- subset(rainfall_deseasonal1, start = 1, end = 3287)
plot(rainfall_train1)
# Test set - 1 year
rainfall_test1 <- subset(rainfall_deseasonal1, start = 3288, end = 3652)
plot(rainfall_test1)
# Fit the model
rainfall_train1_arima <- auto.arima(rainfall_train1, seasonal=TRUE, approximation = FALSE)
# Forecast using the model
rainfall_fcast_test1 <- forecast(rainfall_train1_arima,h=365, seasonal=TRUE, approximation = FALSE)## Warning in forecast.Arima(rainfall_train1_arima, h = 365, seasonal =
## TRUE, : The non-existent seasonal arguments will be ignored.The non-
## existent approximation arguments will be ignored.plot(rainfall_fcast_test1, col="black",xlab="Date", ylab="Precipitation (z-score)")
lines(rainfall_test1,col="red") 
# Evaluating RMSE of forecast
accuracy(rainfall_fcast_test1,rainfall_test1)## ME RMSE MAE MPE MAPE MASE
## Training set -3.918303e-05 0.5315159 0.3574662 86.62683 164.6431 0.6785569
## Test set 6.959079e-02 0.6452657 0.4388499 61.94763 105.3055 0.8330428
## ACF1 Theil's U
## Training set -0.00042333 NA
## Test set 0.45847534 0.9572293Holt-Winters Seasonal Exponential Smoothing
If the data has no trend or seasonal components, simple exponential smoothing can be used for forecasting. If the data has a linear trend component, Holt’s linear method or the damped method can be used. However, if the data is seasonal, as this data is, the Holt-Winters method becomes useful. The Holt-Winters (H-W) method is based on three smoothing equations; one each for level, trend, and seasonality. It is an extension of Holt’s linear method, but there are two H-W methods depending on whether seasonality is modeled in an additive or multiplicative way. The triplet ETS refers to error, trend, and seasonality (in that order) and each of these component is assigned a mode from the selection of possibilities: “additive”, “multiplicative”, or “none”.
If you call forecast directly on a time series object, it will select an ETS model automatically and then return the forecast. Estimation is handled via maximizing the likelihood of the data given to the model, where minimizing the AIC gives the best model for prediction. It works by applying each of 30 smoothing methods that are appropriate to the data, estimating parameters and initial values using MLE, and selecting best method using AIC.
rainfall_ets1 <- stl(rainfall_transformed_ts,s.window="periodic")
plot(forecast(rainfall_ets1, h=365), xlab="Date", ylab="Precipitation (z-score)")
forecast() suggests an exponential smoothing model with error, trend and seasonality terms of “additive”, “none”, and “none” respectively (ETS(A,N,N)).
Compared the other methods tested here, this method reflects more short term dynamics, and produces a more precise forecast.
Testing seasonal exponential smoothing with additive errors
# Training set - 9 years
rainfall_train2 <- subset(rainfall_deseasonal1, start = 1, end = 3287)
# Test set - 1 year
rainfall_test2 <- subset(rainfall_deseasonal1, start = 3288, end = 3652)
# Fit the model
rainfall_train2_ets <- stl(rainfall_train2,s.window="periodic")
# Forecast using the model
rainfall_fcast_test2 <- forecast(rainfall_train2_ets,h=365)
plot(rainfall_fcast_test2, col="black", xlab="Date", ylab="Precipitation (z-score)")
lines(rainfall_test2,col="red") 
# Evaluate RMSE of forecast
accuracy(rainfall_fcast_test2,rainfall_test2)## ME RMSE MAE MPE MAPE
## Training set -0.0004083956 0.5497882 0.3742832 93.631287 207.8522
## Test set 0.2474532717 0.7552154 0.4935506 1.567681 199.6672
## MASE ACF1 Theil's U
## Training set 0.7104796 0.1354041 NA
## Test set 0.9368779 0.4590263 1.351582ARIMA for long seasonality using a Fourier transform as an exogenous covariate
It is a challenge to fit ARIMA or ETS model with data having a long seasonal period, such as 365 for daily data, since seasonal versions of ARIMA and ETS models are designed for shorter periods such as 12 for monthly data or 4 for quarterly data. The problem is that there are m−1 parameters to be estimated for the initial seasonal states where m is the seasonal period. So for large m, the estimation becomes almost impossible.
For such data one can use a Fourier series approach where the seasonal pattern is modelled using Fourier terms with short-term time series dynamics allowed in the error, although in this case the seasonality is assumed to be fixed, i.e. the pattern is not allowed to change over time.
Build ARIMA with fourier term
To do this, we include a Fourier transform in the model as an exogenous covariate. The value of K specifies the number of sine and cosine terms to return for each of the seasonal periods. We wish to capture one annual seasonal cycle, so we execute one Fourier term.
rainfall_arima_fourier <- auto.arima(rainfall_transformed_ts, seasonal=FALSE, approximation = FALSE, xreg=fourier(rainfall_transformed_ts, K=1)) #build model using auto.arima
checkresiduals(rainfall_arima_fourier) #inspect the residuals
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(1,0,2) errors
## Q* = 751.57, df = 724, p-value = 0.2318
##
## Model df: 6. Total lags used: 730plot(forecast(rainfall_arima_fourier, h=365, xreg=fourier(rainfall_transformed_ts, K=1, h=365)), xlab="Date", ylab="Precipitation (z-score)") #plot the resulting forecast
This forecast looks like a pretty good abstraction of the basic annual seasonal cycle projected into the future.
Testing ARIMA for long seasonality using fourier transform as exogenous covariate.
# Training set - 9 years
rainfall_train3 <- subset(rainfall_transformed_ts, start = 1, end = 3287)
# Test set - 1 year
rainfall_test3 <- subset(rainfall_transformed_ts, start = 3288, end = 3652)
# Fit the model
rainfall_train_fourier <- auto.arima(rainfall_train3, approximation = FALSE, seasonal=FALSE, xreg=fourier(rainfall_train3, K=1))
# Forecast using the model
rainfall_fcast_test3 <- forecast(rainfall_train_fourier,h=365, xreg=fourier(rainfall_test3 , K=1))
plot(rainfall_fcast_test3, col="black", xlab="Date", ylab="Precipitation (z-score)")
lines(rainfall_test3,col="red") 
# Evaluate RMSE of forecast
accuracy(rainfall_fcast_test3,rainfall_test3)## ME RMSE MAE MPE MAPE MASE
## Training set -4.525585e-05 0.5644987 0.3833185 60.34252 111.3099 0.7276307
## Test set 6.956225e-02 0.6833967 0.4725633 72.53777 94.7076 0.8970388
## ACF1 Theil's U
## Training set 0.000870831 NA
## Test set 0.431542013 0.8113515Christiano-Fitzgerald band pass filter
The Christiano-Fitzgerald filter is a finite data approximation to the ideal bandpass filter. Several band-pass approximation strategies can be selected in the function cffilter. The default setting of cffilter returns the filtered data yˆt associated with the unrestricted optimal filter assuming no unit root, no drift and an iid filter.
Band pass filters work by combining high pass filters and low pass filters in order to pass frequencies that occur in a defined band or range. High frequency signal components occur when there is relatively small change in value between successive data points; low frequency signal components are occur when there is a relatively large change in value between successive data points. Here we use the filter to remove both high-frequency information and low frequency information in order to isolate two distinct signal components: one corresponding to annual seasonality, and another for long term trend.
Isolate trend-cycle
If the frequency we want to isolate corresponds to the length of the total measurement period, we want to filter out faster frequencies (changes occurring over the course of any period of time smaller than a year) and set the lower pass limit as low as possible. A band pass range of 366-3652 should result in a lower-range band that corresponds to frequencies as high as 1 year long and as low as 10 years long, i.e. the domain of the entire measurement period.
rainfall_cffilter_trend <- cffilter(rainfall$prcp_log,pl=366,pu=3652,root=FALSE,drift=FALSE,
type=c("asymmetric"),
nfix=1,theta=1)
plot(rainfall_cffilter_trend) 
Subtract trend-cycle
We remove the trend signal component before isolating the seasonal signal component.
rainfall_detrended <- (rainfall$prcp_log - rainfall_cffilter_trend$cycle)Isolate seasonal-cycle from de-trended data
If the frequency we want to isolate corresponds to the length of an annual season, we want to filter out faster frequencies (changes occurring over the course of days or weeks) and slower frequencies (gradual change over the course of multiple years). A band pass range of 180-365 should result in a higher-range band that corresponds to frequencies as high as 6 months long and as low as 1 year long, i.e. the domain of a single annual season.
rainfall_cffilter_seas <- cffilter(rainfall_detrended,pl=180,pu=365,root=FALSE,drift=FALSE,
type=c("asymmetric"),
nfix=1,theta=1)
plot(rainfall_cffilter_seas) 
Vizuale trend and seasonal cycle frequencies
plot(rainfall$prcp_log, xlab="Days", ylab="Precipitation")
lines(rainfall_cffilter_trend$cycle, col="red", lwd = 2)
lines(rainfall_cffilter_seas$cycle, col="blue", lwd = 2)
legend(1, 3, legend=c("Trend", "Season"),
col=c("red", "blue"), lty=1:1, cex=1, lwd = 2)
Forecasting with CF Band Pass Method
# convenience dataframe containing trend, seasonal, and date vectors
rainfall_cycles <- data.frame(Date=as.Date(rainfall$water_date),
Seasonal=rainfall_cffilter_seas$cycle,
Trend=rainfall_cffilter_trend$cycle) Inspect trend for stationarity
We do not believe we need the trend-cycle extracted using the low-frequency band pass filter in order to construct the average annual seasonal cycle and forecast it (because we already concluded that the series is trend-stationary and therefore the trend will not influence the forecast), but we can use this different approach to analysis to test our assumption that the trend is stationary.
ggplot(rainfall_cycles, aes(Date, Trend)) +
geom_point() +
scale_x_date(date_labels = "%Y", date_breaks = "1 year", minor_breaks = NULL) +
theme(axis.text.x = element_text(hjust = 0)) +
geom_smooth(alpha = 0.25, size=1, method = "lm", se = FALSE) # Plot the trend-cycle 
adf.test(rainfall_cycles$Trend)##
## Augmented Dickey-Fuller Test
##
## data: rainfall_cycles$Trend
## Dickey-Fuller = -64.457, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationarykpss.test(rainfall_cycles$Trend) #p-values of < 0.05 indicate the trend is stationary##
## KPSS Test for Level Stationarity
##
## data: rainfall_cycles$Trend
## KPSS Level = 0.23692, Truncation lag parameter = 13, p-value = 0.1Our assumption that the trend is stationary is confirmed.
Construct 10-year average annual seasonal cycle
The forecasting technique associated with this method will consists of identifying the 10 exisitng annual-seasonal cycles, regressing them to an average shape, and projecting this as a forecast.
# Create a vector composed of the seasonal cycle for each of 10 year averaged into one
rainfall_seas_avg <- rainfall_cycles %>%
mutate(date_hydro = ymd(paste(
if_else(month(Date) < 10, 2017, 2016),
month(Date), day(Date)))) %>%
group_by(date_hydro) %>%
summarize(multiyearmean = mean(Seasonal))
# Backtransform standardization and normalization
rainfall_seas_backtrans <- rainfall_seas_avg
# (inverse normalization) inverse of the inverse hyperbolic sine function is the hyperbolic sine function
rainfall_seas_backtrans$multiyearmean_sinh <- sinh(rainfall_seas_backtrans$multiyearmean)
# (inverse standardization) backtransforming log using the attributes of original scale() output
rainfall_seas_backtrans$multiyearmean_unscale <- rainfall_seas_backtrans$multiyearmean_sinh * attr(rainfall$prcp_scaled, 'scaled:scale') + attr(rainfall$prcp_scaled, 'scaled:center')
# Plot of backtransformed seasonal avg
ggplot(rainfall_seas_backtrans, aes(date_hydro, multiyearmean_unscale)) +
geom_point() +
scale_x_date(date_labels = "%b", date_breaks = "1 month", minor_breaks = NULL) +
theme(axis.text.x = element_text(hjust = 0))
Testing the bandpass-derived forecast
# Convenience object containing only time and precipitation variables
rainfall_observations <- rainfall %>% ungroup() %>% select(water_date, PRCP)
# Create training and testing sets
rainfall_train4 <- rainfall_observations[c(1:3288),c(1:2)]
rainfall_test4 <- rainfall_observations[c(3289:3653),c(1:2)]
# Create seasonal average
rainfall_bandpass_forecast <- rainfall_seas_backtrans %>% select(date_hydro, multiyearmean_unscale) %>% 'colnames<-'(c("water_date", "PRCP"))
# Visualize forecast compared to observations
plot(rainfall_train4, col="black")
lines(rainfall_test4,col="red")
lines(rainfall_bandpass_forecast,col="blue")
# Test forecast accuracy
rainfall_bandapass_fcast_ts <- ts(rainfall_bandpass_forecast$PRCP, start=c(2016, 10), end=c(2017, 9), frequency = 365) #forecast
rainfall_observed_ts <- ts(rainfall_test4$PRCP, start=c(2016, 10), end=c(2017, 9), frequency = 365) #observed data from forecast period
accuracy(rainfall_bandapass_fcast_ts,rainfall_observed_ts)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.04391553 0.4515508 0.2680407 -Inf Inf 0.3722014 0Forecast Accuracy Comparison
We now compare the accuracy of each forecast method. The accuracy() function of the forecast package will return accuracy evaluations for both in-sample and out-of-sample data if both are provided.
## Seasonal ARIMA
arima1_acc <- accuracy(rainfall_fcast_test1,rainfall_test1)
arima1_acc## ME RMSE MAE MPE MAPE MASE
## Training set -3.918303e-05 0.5315159 0.3574662 86.62683 164.6431 0.6785569
## Test set 6.959079e-02 0.6452657 0.4388499 61.94763 105.3055 0.8330428
## ACF1 Theil's U
## Training set -0.00042333 NA
## Test set 0.45847534 0.9572293## Holt-Winters Exponential Smoothing
ets_acc <- accuracy(rainfall_fcast_test2,rainfall_test2)
ets_acc## ME RMSE MAE MPE MAPE
## Training set -0.0004083956 0.5497882 0.3742832 93.631287 207.8522
## Test set 0.2474532717 0.7552154 0.4935506 1.567681 199.6672
## MASE ACF1 Theil's U
## Training set 0.7104796 0.1354041 NA
## Test set 0.9368779 0.4590263 1.351582## Non-Seasonal ARIMA with Fourier Term
arima2_acc <- accuracy(rainfall_fcast_test3,rainfall_test3)
arima2_acc## ME RMSE MAE MPE MAPE MASE
## Training set -4.525585e-05 0.5644987 0.3833185 60.34252 111.3099 0.7276307
## Test set 6.956225e-02 0.6833967 0.4725633 72.53777 94.7076 0.8970388
## ACF1 Theil's U
## Training set 0.000870831 NA
## Test set 0.431542013 0.8113515## Christiano-Fitzgerlad Bandpass Filter
cffilter_acc <- accuracy(rainfall_bandapass_fcast_ts,rainfall_observed_ts)
cffilter_acc## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.04391553 0.4515508 0.2680407 -Inf Inf 0.3722014 0Plot the accuracy results
From the plot we can see that the Holt-Winters Seasonal Exponential Smoothing method outperformed all of the other methods on both the in-sample and out-of-sample data sets, except for the ARIMA with Fourier term method, which outperformed it on the out-of-sample data set only. H-W exponential smoothing also admits of the largest performance discrepancy between in-sample and out-of-sample forecast performance, of the three pairs recorded.
#convert a`ccuracy()` results to dataframes
acc1 = data.frame(arima1_acc)
acc2 = data.frame(ets_acc)
acc3 = data.frame(arima2_acc)
acc4 = data.frame(cffilter_acc)
acc4$MASE <- NA
acc4 <- rbind(acc4, "Training set" = c(NA,NA,NA,NA,NA,NA,NA,NA))
#combine dataframes into a list
lsAcc = list(acc1, acc2, acc3, acc4)
dfAcc = do.call(what = rbind, args = lsAcc)
dfAcc$rowPair = ceiling(1:8 * 0.5)
#plot training set and testing set RMSE for each method
plot(0, cex=0, xlim=c(1,4), ylim=c(0,1), xlab="Methods", ylab="RMSE Score")
points(x = dfAcc$rowPair[grepl(pattern = "Training", rownames(dfAcc))],
y = dfAcc$RMSE[grepl(pattern = "Training", rownames(dfAcc))],
pch = 2, col="blue")
points(x = dfAcc$rowPair[grepl(pattern = "Test", rownames(dfAcc))],
y = dfAcc$RMSE[grepl(pattern = "Test", rownames(dfAcc))],
pch = 3, col="red")
Conclusion
This experiment is concluded by asserting that, of the four methods tested here, based on the data utilized and the assumptions held in the preceding arguments, the Holt-Winters Seasonal Exponential Smoothing technique produced the overall most accurate and appropriate 1-year-out forecast of precipitation levels.
Areas for Improvement
The forecast produced by HW Seasonal Exponential Smoothing did achieve the highest RMSE score. Not only that, the shape of its forecast is the most morphologically realistic of the tested options. However, intuition tells us that while it produces a very precise forecast, this degree of precision will not necessarily yield the most accurate measurements at a given point in time. A simpler regression might yield less precise forecast while achieving a more accurate prediction for any given measurement.
The HW method may be yielding an over-fitted model, which may be what is causing the inappropriately precise forecast. Over-fitting can be a problem for univariate time series prediction because there are no additional dimensions of data to balance the patterns detected in the single variable and inform the forecast.
K-fold Cross Validation should be explored as a method of testing that is superior to the single holdout validation set method. K-fold cross-validation could estimate the model’s predictive performance by repeating the holdout experiment again for each year of data, or even within windows of time that do not correspond with single annual seasons. This may result in a less favorable evaluation of the HW method’s possibly overfitted prediction.
Other measurements of goodness-of-fit besides RMSE, such as mean absolute scaled error (MASE) could also be evaluated and may yield meaningfully different results.
Analysing a larger data set may also improve the accuracy of these conclusions. For example, over a longer time scale, non-stationarity may become evident.
It is also worth noting that the characteristic patterns of the rainfall recorded in this data set are not represented across all climates and areas and that these techniques may perform differently in different contexts. It would be instructive to compare the results on these same analysis and forecasting methods across multiple, morphologically dissimilar precipitation data sets.