Redes neurais de séries temporais no R

Uma rede neural artificial é um método de previsão baseado na modelagem matemática do cérebro: os neurônios são organizados em camadas de entrada (preditores) e uma camada de saída (previsão), podendo haver camadas intermediárias/ocultas. O modo mais simples seria uma regressão linear, cujos regressores (variáveis independentes) são multiplicados por coeficientes (pesos), que são reajustados a cada loop de modo a minimizar o erro. Ao adicionar camadas ocultas, a rede torna-se não-linear. Veja mais sobre o que são redes neurais artificias clicando no link.

Os valores atrasados de uma série temporal (inclusive com sazonalidade) podem ser usados como entradas para uma rede neural. Considerando-se uma rede unidirecional (feedfoward) com uma camada oculta, utiliza-se a notação NNAR(p,P,k)[m] para indicar que o atraso (lag) é de p observações (por exemplo, para p=9, são as 9 últimas observações), com ordens de sazonalidade P e m (por exemplo, para P=1 e m=12, é considerado o valor de 12 amostras/meses atrás) e k nós (neurônios) na camada oculta.

A função nnetar() do pacote forecast do R ajusta um modelo do modo apresentado acima. Para uma série não-sazonal, o valor padrão de p é definido pelo número ótimo de uma AR(p), conforme dado pelo AIC; para uma série sazonal, o valor padrão de P é 1, p é escolhido a partir do modelo linear ótimo ajustado aos dados sazonais e o valor inteiro arredondado da conta k=(p+P+1)/2.

Se rodar o cálculo da rede neural (função nnetar) mais de uma vez, os resultados saem diferentes (já que a inicialização dos pesos é dada com valores aleatórios), mas todos estão dentro do intervalo de confiança. É possível definir a semente randômica usando “set.seed(12345)” logo antes de chamar a nnetar. No entanto, salvando os resultados e fazendo uma média, o valor resultante da previsão não varia muito se rodar o script várias vezes.

Assim, em vez de usar diretamente a função forecast() aplicada à rede neural, pode-se usar forecast.nnetar() para controlar quantas simulações serão feitas (o padrão é 1000, mas pode ser especificado outro valor usando o parâmetro “npaths=”). Outros parâmetros podem ser especificados, como “repeats”, que dá o número de redes calculadas com diferentes pesos iniciais aleatórios e que serão calculados como média quando produzem previsões; “decay” é a taxa de aprendizagem e “maxit” é o máximo de interações permitidas.

O script a seguir pode ser copiado para um arquivo e salvo como “rede_neural.R”, sendo executado diretamente no terminal (./rede_neural.R) após ter sua permissão de execução incluída (chmod +x rede_neural.R). Ele abre os arquivos com a variável a ser prevista e os regressores (passados e futuros), além de um período usado para testes (os valores são conhecidos, mas não fazem parte do grupo de treinamento da rede). Foram construídas séries temporais para entrada no modelo e na previsão para 12 meses (h=12) – intervalos de confiança podem ser calculados incluindo “PI=TRUE”. Pode-se utilizar uma transformação Box-Cox (lambda) para garantir que os resíduos sejam aproximadamente homoscedásticos.

#!/usr/bin/Rscript
# Script para calcular previsão de série temporal usando redes neurais

# Definindo bibliotecas
end_libs="~/Rpacks"
suppressPackageStartupMessages(require(timeDate,lib=end_libs))
suppressPackageStartupMessages(require(zoo,lib=end_libs))
suppressPackageStartupMessages(require(forecast,lib=end_libs))
options(warn = 1) # If 'warn' is one, warnings are printed as they occur

# Obter séries de dados
var_obs = read.csv('dataset2.csv', header = TRUE)
reg_obs = read.csv('xregobs2.csv', header = TRUE)
reg_prev = read.csv('xregprev2.csv', header = TRUE)
resp = read.csv('respostas2.csv', header = TRUE)

# Construir séries temporais
temp = strsplit(as.character(var_obs$data[1]), split="-")
start_year_hist = as.numeric(temp[[1]][1])
start_month_hist = as.numeric(temp[[1]][2])
varobs_ts = ts(var_obs[,'valor'], frequency = 12, start = c(start_year_hist, start_month_hist))
regobs_ts = ts(reg_obs[,'valor'], frequency = 12, start = c(start_year_hist, start_month_hist))
temp = strsplit(as.character(reg_prev$data[1]), split="-")
start_year_prev = as.numeric(temp[[1]][1])
start_month_prev = as.numeric(temp[[1]][2])
regprev_ts = ts(reg_prev[,'valor'], frequency = 12, start = c(start_year_prev, start_month_prev))
respostas = ts(resp[,'valor'], frequency = 12, start = c(start_year_prev, start_month_prev))

# Criar modelo
lambda = BoxCox.lambda(varobs_ts)
fit = nnetar(varobs_ts, xreg = regobs_ts, decay=0.001, repeats=1000, lambda = lambda)
# Fazer previsão
previsao = forecast.nnetar(fit, xreg = regprev_ts, h = 12, lambda = lambda)

# Imprimir/plotar/gravar resultados
print(round(previsao$mean,digits=0))
print(round(respostas,digits=0))
print(round(previsao$mean,digits=0)-round(respostas,digits=0))
plot(previsao)

A saída imprime (na sequência) o valor previsto, o esperado e a diferença entre eles. Note que a série ainda precisa receber mais tratamentos. através da aplicação de outras técnicas, para gerar valores melhores.

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2016                      79  90  88 110 107 130 135
2017 139 157 138 121 112                            
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2016                     121  92 104  94  85  89  79
2017  69  70  69  73 101                            
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2016                     -42  -2 -16  16  22  41  56
2017  70  87  69  48  11     

Outra opção é o multilayer perceptrons (MLP).

Fontes

• Rob J Hyndman – Prediction intervals for NNETAR models
• Rob J Hyndman & George Athana­sopou­los – Neural network models

Posts relacionados:

Estatística multivariada Tendência e sazonalidade Modelo ARIMA e previsão no R Correlação