Nonparametric Risk Bounds for Time-Series Forecasting

Abstract: We derive generalization error bounds for traditional time- series forecasting models. Our results hold for many standard forecasting tools including autoregressive models, moving average models, and, more generally, linear state-space models. These non-asymptotic bounds need only weak assumptions on the data-generating process, yet allow forecasters to select among competing models and to guarantee, with high probability, that their chosen model will perform well. We motivate our techniques with and apply them to standard economic and financial forecasting tools—a GARCH model for predicting equity volatility and a dynamic stochastic general equilibrium model (DSGE), the standard tool in macroeconomic forecasting. We demonstrate in particular how our techniques can aid forecasters and policy makers in choosing models which behave well under uncertainty and mis-specification.
Conclusion: This paper demonstrates how to control the generalization error of common time-series forecasting models, especially those used in economics and engineering—ARMA models, vector autoregressions (Bayesian or otherwise), linearized dynamic stochastic general equilibrium models, and linear state-space models. We derive upper bounds on the risk, which hold with high probability while requiring only weak assumptions on the data-generating process. These bounds are finite sample in nature, unlike standard model selection penalties such as AIC or BIC. Furthermore, they do not suffer the biases inherent in other risk estimation techniques such as the pseudo-cross validation approach often used in the economic forecasting literature.
While we have stated these results in terms of standard economic forecasting models, they have very wide applicability. Theorem 12 applies to any forecasting procedure with fixed memory length, linear or non-linear. Theorem 17 applies only to methods whose forecasts are linear in the observations, but a similar result for nonlinear methods would just need to ensure that the dependence of the forecast on the past decays in some suitable way. Rather than deriving bounds theoretically, one could attempt to estimate bounds on the risk. While cross-validation is tricky (Racine, 2000), nonparametric bootstrap procedures may do better. A fully nonparametric version is possible, using the circular bootstrap (reviewed in Lahiri, 1999). Bootstrapping lengthy out-of-sample sequences for testing fitted model predictions yields intuitively sensible estimates of Rn(f), but there is currently no theory about the coverage level. Also, while models like VARs can be fit quickly to simulated data, general state-space models, let alone DSGEs, require large amounts of computational power, which is an obstacle to any resampling method.
Nonparametric Risk Bounds for Time-Series Forecasting

Previsão do número de homicídios no Rio de Janeiro (Sim, eu voltei a postar!!!)

Um pequeno disclaimer: Depois de muito relutar, estou voltando com as atividades do blog. Isso significa de maneira prática, que apesar do clipping ser o grande driver dos acessos, sinto que há uma lacuna em relação a posts técnicos aqui na blogosfera brasileira.

Dessa forma, vou fazer mais posts técnicos com assuntos ligados à Machine Learning, Data Mining e afins.

Disclaimer feito, bora pro post original.

*******

Fala pessoal, tudo bem?

Um dos principais problemas em utilizar processadores numéricos como o Excel ou mesmo o LibreOffice é que se perde um pouco a flexibilidade de parametrizar os modelos de previsão.

Além disso, há um problema grave nesses processadores de reprodutibilidade, dado que pode haver inúmeras versões, planilhas voando pra todo o lado via e-mail, e por aí vai.

Nesse post eu vou utilizar dados reais da quantidade de homicídios que houveram no Rio de Janeiro de 1991 até 2011.

Os dados foram extraídos e consolidados pelo CESeC – Centro de Estudos de Segurança e Cidadania e o link está no final do script.

Antes de mais nada, vamos carregar o package Forecast:

# Load library
if(!require(forecast)){
 install.packages("forecast")
 library(forecast)
}

Depois disso, vamos realizar a carga dos dados. Nesse caso eu fiz uma transformação na planilha original e passei os dados para csv.

# Dataset
rj_murder <- read.csv("https://raw.githubusercontent.com/fclesio/learning-space/master/Datasets/rj-homicidios-1991-2011.csv")

Vamos dar uma olhada nos dados para nos certificar que tudo está OK.

# Check data
head(rj_murder)

Com os dados OK, para incluir um conjunto de dados no package de Forecast e ARIMA, temos que antes de mais nada, pegar todos os datapoints e construir uma série através de um objeto de timeseries como no bloco de código abaixo:

# Transform in time series dataset all metrics inside the dataset 
rj_murder_total <- ts(rj_murder$total, start=c(1991, 1), end=c(2011, 12), frequency=12)
rj_murder_capital <- ts(rj_murder$capital, start=c(1991, 1), end=c(2011, 12), frequency=12)
rj_murder_baixada <- ts(rj_murder$baixada, start=c(1991, 1), end=c(2011, 12), frequency=12)
rj_murder_interior <- ts(rj_murder$interior, start=c(1991, 1), end=c(2011, 12), frequency=12)

Para entendermos o que aconteceu, vamos destrinchar erssa função TS:

ts() função para a construção de série temporal

rj_murder$total conjunto de datapoints do dataframe da coluna total

start=c(1991, 1) início da série é janeiro de 1991

end=c(2011, 12) o fim da série é dezebro de 2011

frequency=12 como estamos falando de dados mensais, a frequencia da série é de 12 datapoints por ano

Agora que entendemos o que fizemos, vamos dar uma olhada nos nossos objetos ts() e ver se os mesmos estão corretos.

#Plot series
par(mfrow=c(2,2))

plot(rj_murder_total
 ,main = "Total de Homicidios no RJ de 1991-2011"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 )


plot(rj_murder_capital
 ,main = "Total de Homicidios no RJ de 1991-2011 (Capital)"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
)


plot(rj_murder_baixada
 ,main = "Total de Homicidios no RJ de 1991-2011 (Baixada)"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
)


plot(rj_murder_interior
 ,main = "Total de Homicidios no RJ de 1991-2011 (Interior)"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
)

 

Homicidios RJ

 

Agora que vimos que as nossas séries estão corretas, vamos criar alguns objetos usando a função auto.arima().

# Fit ARIMA models
fit_rj_murder_total <- auto.arima(rj_murder_total)
fit_rj_murder_capital <- auto.arima(rj_murder_capital) 
fit_rj_murder_baixada <- auto.arima(rj_murder_baixada)
fit_rj_murder_interior <- auto.arima(rj_murder_interior)

Essa função retorna o melhor modelo ARIMA de acordo com os valores de AIC, AICc e BIC.

Mas deixando um pouco a teoria de lado, vamos ver como ficaram os ajustes dos modelos ARIMA de acordo com cada uma das séries.

#Plot ARIMA Models
par(mfrow=c(2,2))
plot(forecast(fit_rj_murder_total,h=12)
 ,main = "Total de Homicidios no RJ de 1991-2011 \n Previsão usando ARIMA"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5)

plot(forecast(fit_rj_murder_capital,h=12)
 ,main = "Total de Homicidios no RJ de 91-2011 (Capital) \n Previsão usando ARIMA"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5)

plot(forecast(fit_rj_murder_baixada,h=12)
 ,main = "Total de Homicidios no RJ de 91-2011 (Baixada) \n Previsão usando ARIMA"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5)

plot(forecast(fit_rj_murder_interior,h=12)
 ,main = "Total de Homicidios no RJ de 91-2011 (Interior) \n Previsão usando ARIMA"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5)

 

ARIMA

 

Um bom resultado, mas agora vamos ver o que conseguimos fazer usando o package forecast default, com o horizonte de previsão de 12 meses:

#Forecasting using confidence intervals in 80%, 95% and 99% with 12 months ahead
ahead_rj_murder_total <- forecast(rj_murder_total, level=c(80, 95, 99), h=12)
ahead_rj_murder_capital <- forecast(rj_murder_capital, level=c(80, 95, 99), h=12)
ahead_rj_murder_interior <- forecast(rj_murder_interior, level=c(80, 95, 99), h=12)
ahead_rj_murder_baixada <- forecast(rj_murder_baixada, level=c(80, 95, 99), h=12)

Modelos ajustados, vamos plotar os resultados:

#Plot forecasting
par(mfrow=c(2,2))
plot(ahead_rj_murder_total
 ,main = "Total de Homicidios no RJ de 91-2011 (Total) \n Previsão usando Forecast package"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 ,shadecols="oldstyle")

plot(ahead_rj_murder_capital
 ,main = "Total de Homicidios no RJ de 91-2011 (Capital) \n Previsão usando Forecast package"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 ,shadecols="oldstyle")

plot(ahead_rj_murder_baixada
 ,main = "Total de Homicidios no RJ de 91-2011 (Baixada) \n Previsão usando Forecast package"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 ,shadecols="oldstyle")

plot(ahead_rj_murder_interior
 ,main = "Total de Homicidios no RJ de 91-2011 (Interior) \n Previsão usando Forecast package"
 ,xlab = "Ano"
 ,ylab = "Qtde Homicidios"
 ,type = "l"
 ,lwd = 2.5
 ,shadecols="oldstyle")

Forecast Package Predictions

Podemos ver que tanto usando o modelo ARIMA quanto com o uso do package forecast vimos que haverá uma queda no número de homicídios. Com isso podem ser elaboradas (usando mais dados) políticas públicas de segurança.

Refs:

[1] – http://www.inside-r.org/packages/cran/forecast/docs/auto.arima

[2] – http://www.ucamcesec.com.br/wordpress/wp-content/uploads/2011/04/RJ_2011_evolucao.xls

[3] – http://www.statmethods.net/advstats/timeseries.html

 

Previsão do número de homicídios no Rio de Janeiro (Sim, eu voltei a postar!!!)