Time Series Forecasting
This example shows time series forecasting of Euro-AUD exchange rates with the with the ARIMA and STL models. The data used are historical currency exchange rates from January 1999 to June 2014 provided by the European Central Bank.
This example was produced with R Markdown. The Rmd and R source code files are provided at the bottom of this page.
1. Downloading data from European Central Bank
Download data from the European Central Bank at http://www.ecb.europa.eu/stats/exchange/eurofxref/html/index.en.html.
url <- "http://www.ecb.europa.eu/stats/eurofxref/eurofxref-hist.zip" download.file(url, "eurofxref-hist.zip")
2. Checking data
rates <- read.csv(unz("eurofxref-hist.zip", "eurofxref-hist.csv"),
header = T)rates[1:2, ]
## Date USD JPY BGN CYP CZK DKK EEK GBP HUF LTL LVL ## 1 2014-07-01 1.369 139.0 1.9558 N/A 27.43 7.456 N/A 0.7981 310.4 3.453 N/A ## 2 2014-06-30 1.366 138.4 1.9558 N/A 27.45 7.456 N/A 0.8015 309.3 3.453 N/A ## MTL PLN ROL RON SEK SIT SKK CHF ISK NOK HRK RUB TRL ## 1 N/A 4.158 N/A 4.3881 9.160 N/A N/A 1.214 N/A 8.438 7.58 46.895 N/A ## 2 N/A 4.157 N/A 4.383 9.176 N/A N/A 1.216 N/A 8.403 7.576 46.3779 N/A ## TRY AUD BRL CAD CNY HKD IDR INR KRW MXN ## 1 2.9066 1.448 3.0349 1.459 8.4883 10.61 16251.94 82.2307 1385 17.7759 ## 2 2.8969 1.454 3.0002 1.459 8.4722 10.59 16248.15 82.2023 1382 17.7124 ## MYR NZD PHP SGD THB ZAR ILS X ## 1 4.3893 1.562 59.764 1.705 44.367 14.58 4.692 NA ## 2 4.3856 1.563 59.652 1.705 44.323 14.46 4.696 NA
str(rates$Date)
## Factor w/ 3968 levels "1999-01-04","1999-01-05",..: 3968 3967 3966 3965 3964 3963 3962 3961 3960 3959 ...
## convert into date formatrates$Date <- as.Date(rates$Date, "%Y-%m-%d")str(rates$Date)
## Date[1:3968], format: "2014-07-01" "2014-06-30" "2014-06-27" "2014-06-26" ...
range(rates$Date)
## [1] "1999-01-04" "2014-07-01"
rates <- rates[order(rates$Date), ]
## plot time seriesplot(rates$Date, rates$AUD, type = "l")
3. Forecasting with ARIMA
The code below shows that there are no data for weekends or public holidays.
head(rates$Date, 20)
## [1] "1999-01-04" "1999-01-05" "1999-01-06" "1999-01-07" "1999-01-08" ## [6] "1999-01-11" "1999-01-12" "1999-01-13" "1999-01-14" "1999-01-15" ## [11] "1999-01-18" "1999-01-19" "1999-01-20" "1999-01-21" "1999-01-22" ## [16] "1999-01-25" "1999-01-26" "1999-01-27" "1999-01-28" "1999-01-29"
years <- format(rates$Date, "%Y")tab <- table(years)tab
## years ## 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 ## 259 255 254 255 255 259 257 255 255 256 256 258 257 256 255 ## 2014 ## 126
## number of days per year after removing 2014mean(tab[1:(length(tab) - 1)])
## [1] 256.1
Based on above result, there are about 256 values per year, so the windows size is set to 256 in time series analysis in section 5. Another way is to fill weekends and public holidays with values in the previous populated days.
source("forecast.R") ## see code file in section 5result.arima <- forecastArima(rates, n.ahead = 90)
source("plotForecastResult.R") ## see code file in section 5plotForecastResult(result.arima, title = "Exchange rate forecasting with ARIMA")
4. Forecasting with STL
result.stl <- forecastStl(rates, n.ahead = 90)
plotForecastResult(result.stl, title = "Exchange rate forecasting with STL")
## exchange rate in 2014result <- subset(result.stl, date >= "2014-01-01")plotForecastResult(result, title = "Exchange rate forecasting with STL (2014)")
5. Functions
Below are two source files used in section 3 and 4.
It provides functions for forecasting with ARIMA and STL.
>
library(forecast) > forecastStl <- function(x, n.ahead = 30) { + myTs <- ts(x$AUD, start = 1, frequency = 256) + fit.stl <- stl(myTs, s.window = 256) + sts <- fit.stl$time.series + trend <- sts[, "trend"] + fore <- forecast(fit.stl, h = n.ahead, level = 95) + plot(fore) + pred <- fore$mean + upper <- fore$upper + lower <- fore$lower + output <- data.frame(actual = c(x$AUD, rep(NA, n.ahead)), + trend = c(trend, rep(NA, n.ahead)), pred = c(rep(NA, + nrow(x)), pred), lower = c(rep(NA, nrow(x)), lower), + upper = c(rep(NA, nrow(x)), upper), date = c(x$Date, + max(x$Date) + (1:n.ahead))) + return(output) + } > forecastArima <- function(x, n.ahead = 30) { + myTs <- ts(x$AUD, start = 1, frequency = 256) + fit.arima <- arima(myTs, order = c(0, 0, 1)) + fore <- forecast(fit.arima, h = n.ahead) + plot(fore) + upper <- fore$upper[, "95%"] + lower <- fore$lower[, "95%"] + trend <- as.numeric(fore$fitted) + pred <- as.numeric(fore$mean) + output <- data.frame(actual = c(x$AUD, rep(NA, n.ahead)), + trend = c(trend, rep(NA, n.ahead)), pred = c(rep(NA, + nrow(x)), pred), lower = c(rep(NA, nrow(x)), lower), + upper = c(rep(NA, nrow(x)), upper), date = c(x$Date, + max(x$Date) + (1:n.ahead))) + return(output) + }
It provides a function for ploting time series forecasting result, incl. trend, forecast and bounds.
> plotForecastResult <- function(x, title = NULL) { + x <- x[order(x$date), ] + max.val <- max(c(x$actual, x$upper), na.rm = T) + min.val <- min(c(x$actual, x$lower), na.rm = T) + plot(x$date, x$actual, type = "l", col = "grey", main = title, + xlab = "Time", ylab = "Exchange Rate", xlim = range(x$date), + ylim = c(min.val, max.val)) + grid() + lines(x$date, x$trend, col = "yellowgreen") + lines(x$date, x$pred, col = "green") + lines(x$date, x$lower, col = "blue") + lines(x$date, x$upper, col = "blue") + legend("bottomleft", col = c("grey", "yellowgreen", "green", + "blue"), lty = 1, c("Actual", "Trend", "Forecast", "Lower/Upper Bound")) + }