Tag

modelo sarima Archives - Análise Macro

Um modelo univariado para Homicídios por Intervenção Policial no RJ

By | Clube AM

Ontem, o Jornal Nacional apresentou uma matéria sobre o número recorde de homicídios por intervenção policial no Rio de Janeiro nesse ano, dando como causa a nova política de segurança pública do governador Wilson Witzel. Fiquei curioso e fui dar uma olhada na ótima base de dados do Instituto de Segurança Pública (ISP). Peguei o arquivo Séries históricas do estado por mês desde 1991 (números absolutos) e importei para o R com o código abaixo, já carregando os pacotes que iria utilizar na análise.


library(readr)
library(forecast)
library(ggplot2)
library(TStools)
library(scales)

data = read_csv2('DOMensalEstadoDesde1991.csv')

 

Os dados de homicídios por intervenção policial estão disponíveis desde janeiro de 1998, apenas. Com isso, fazemos uma amostra dos nossos dados com o código a seguir.


data = ts(data, start=c(1991,01), freq=12)
hip = window(data[,6], start=c(1998,01))

E um gráfico rápido para visualizar os nossos dados.


df = data.frame(time=as.Date(time(hip)),
hip=hip, cmav=cmav(hip,12))

ggplot(df, aes(x=time))+
geom_line(aes(y=hip, colour='HIP'))+
geom_line(aes(y=cmav, colour='Média Móvel'), size=.8)+
scale_colour_manual('',
values=c('HIP'='black',
'Média Móvel'='red'))+
theme(legend.position = c(.15,.8))+
theme(axis.text.x = element_text(angle=45, hjust=1))+
scale_x_date(breaks = date_breaks("1 year"),
labels = date_format("%Y"))+
labs(x='', y='quantidade de pessoas',
title='Homícidios por intervenção policial no Rio de Janeiro',
caption='Fonte: analisemacro.com.br com dados do ISP/RJ')

Observe que acrescentei uma média móvel de 12 meses nos dados, de modo a suavizar a série. O que se observa é de fato um aumento preocupante na série desde ao menos 2013. Ou seja, o aumento é bem anterior à gestão do atual governador. A matéria do JN está, portanto, equivocada ao apontar como causa a nova política de segurança.

Como não sou especialista no assunto, não posso opinar sobre as causas para esse aumento. Minha função aqui é apenas ilustrar como tratar os dados. Isso dito, resolvi fazer uma projeção usando um modelo univariado simples, isto é, um modelo ARIMA, usando o pacote forecast. O gráfico abaixo ilustra o ajuste do modelo.

E abaixo um gráfico com a projeção para os próximos 12 meses.

O código completo do exercício estará disponível na Edição 66 do Clube do Código.

Avaliando a acurácia de uma previsão com o R

By | Hackeando o R

Quem trabalha com previsões quantitativas sabe que uma parte importante de qualquer projeto é a avaliação das previsões geradas pelo(s) modelo(s). É dessa forma que podemos verificar se estamos no caminho correto. Para ilustrar como isso pode ser feito no R, vou construir nesse post um modelo univariado simples para a taxa de desemprego medida pela PNAD Contínua e depois efetuar o processo de avaliação das previsões geradas. Para começar, vamos carregar alguns pacotes...


library(forecast)
library(ggplot2)
library(sidrar)
library(xtable)

Vamos pegar os dados que precisamos do SIDRA com o pacote sidrar...


# Dados Brutos
table = get_sidra(api='/t/6318/n1/all/v/1641/p/all/c629/all')
# Pegar a PEA
pea = table$Valor[table$`Condição em relação à força de trabalho e condição de ocupação (Código)`==32386]
# Pegar a População Desocupada
desocupada = table$Valor[table$`Condição em relação à força de trabalho e condição de ocupação (Código)`==32446]
# Criar Desemprego
desemprego = ts(desocupada/pea*100, start=c(2012,03), freq=12)

Com a função ggtsdisplay do pacote forecast podemos visualizar a nossa série e as funções de autocorrelação como abaixo.

Para a construção de um modelo univariado, teríamos que (1) verificar se a nossa série é estacionária e (2) tentar identificar os coeficientes AR e MA através de funções de autocorrelação. Vamos aqui, entretanto, utilizar um algoritmo de modo a automatizar esse processo através de critérios de informação com a função auto.arima:


sarima = auto.arima(desemprego, max.p=2, max.q=4, max.P = 2, max.Q=2)

Aplicada a função auto.arima sobre a nossa série, temos um modelo SARIMA(1,2,0)(1,1,0)[12]. Com base nesse modelo, podemos agora gerar previsões e, consequentemente, avaliá-las. É comum aqui, por suposto, dividir a nossa amostra em duas subamostras: uma de treino, onde rodamos o nosso modelo e outra de teste, onde são feitas as previsões. Nessa amostra de teste é onde ocorre a comparação com as observações efetivas da nossa série. Como regra de bolso, é comum destinar 70% para a subamostra de treino e o restante para a amostra de teste. Como nossa amostra é, entretanto, curta, vou aqui reservar apenas as 6 últimas observações para o conjunto de teste, de modo a ilustrar apenas o código.


training = window(desemprego, end=end(desemprego)-c(0,6))
test = window(desemprego, start=end(desemprego)-c(0,5))
acuracia = Arima(training, order=c(1,2,0), seasonal = c(1,1,0))
acuraciaf = forecast(acuracia, h=length(test), level=95)

Com o código acima, criamos as nossas subamostras, rodamos o modelo na subamostra de treino e gerar a previsão com base no tamanho do vetor de teste. Por fim, com a função accuracy do pacote forecast podemos avaliar essas previsões.


acc = accuracy(acuraciaf$mean, test)

E aí está a tabelinha que queríamos...

ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set -0.19 0.23 0.20 -1.51 1.57 0.37 0.87

O erro médio do nosso modelinho é de 0,19 negativos e o RMSE é de 0,23. Para maiores detalhes sobre essas métricas, veja esse post aqui.

Curtiu o tema? Nós exploramos isso e muito mais no nosso Curso de Construção de Cenários e Previsões usando o R, voltado exclusivamente para a construção de previsões quantitativas no R. Dê uma olhada lá e se inscreva na próxima turma!

Receba diretamente em seu e-mail gratuitamente nossas promoções especiais
e conteúdos exclusivos sobre Análise de Dados!

Assinar Gratuitamente