Business forecasting: Dynamic regression [modelando o PIB de FIJI]
Pedro Costa Ferreira
Vamos conhecer um pouco sobre Fiji?!
A República das Fiji ou das Fiji, por vezes escrito Fidji, em fijiano: Matanitu ko Viti; em inglês: Republic of Fiji, pronunciado, híndi fijiano: फ़िजी गणराज्य Fijī Gaṇarājya), anteriormente conhecida como República das Ilhas Fiji, é um país insular da Oceania, composto por 332 ilhas no Oceano Pacífico.
Fiji faz fronteira marítima com Tuvalu e com o território francês de Wallis e Futuna a norte, com Tonga a leste, e com o território francês da Nova Caledónia, com Vanuatu e com as Ilhas Salomão a oeste. A sul, a terra mais próxima corresponde às ilhas neozelandesas de Kermadec, mas estão bastante afastadas.[wikipidia]
Veja algumas informações interessantes sobre a economia:
Com solos férteis que permitem o cultivo da cana-de-açúcar, este é o principal produto de exportação, além do coco e gengibre.
A economia de Fiji é uma das mais desenvolvidas do Oceano Pacífico. A exportação de açúcar, as remessas de dinheiro feitas por fijianos residentes no exterior, e o turismo (de 300 a 400 mil turistas por ano visitam o país) são as principais fontes de receitas do país.
O açúcar fijiano tem como destino principal a União Europeia. Em 2007 a ajuda econômica da UE ao país foi suspensa, enquanto o governo que assumiu o poder após o golpe de estado não estabelece um cronograma para novas eleições.
Primeiros passos
1 - Setar a pasta em que estamos trabalhando;
2 - ler os dados;
3 - "Comunicar" ao R que temos uma série de tempo;
4 - Análise exploratória.
1 - Setar a pasta em que estamos trabalhando;
setwd("sua pasta de trabalho")
2 - ler os dados;
require(xlsx)## Loading required package: xlsx## Loading required package: rJava## Loading required package: xlsxjarsFiji_dataset <- read.xlsx("Fiji.xlsx", sheetName="Fiji",as.data.frame = TRUE, header = TRUE)
#View(Fiji_dataset)
names(Fiji_dataset)## [1] "Date" "RGDP_FIJI"
## [3] "G_REAL_GROWTH_FIJI" "RGDP_AUS"
## [5] "G_REAL_GROWTH_AUS" "dummy_crise"
## [7] "REER_FIJI" "SUGAROUTPUT_FIJI"
## [9] "G_SUGAR_OUTPUT" "SUGAROUTPUT_FIJI_MISSING"
## [11] "NA."attach(Fiji_dataset)
3 - “Comunicar” ao R que temos uma série de tempo;
RGDP_FIJI <- ts(RGDP_FIJI,start = 1970,end = 2013,frequency = 1)
RGDP_AUS <- ts(RGDP_AUS,start = 1970,end = 2013,frequency = 1)
SUGAROUTPUT_FIJI <- ts(SUGAROUTPUT_FIJI,start = 1970,end = 2013,frequency = 1)
G_REAL_GROWTH_FIJI <- ts(G_REAL_GROWTH_FIJI,start = 1970,end = 2013,frequency = 1)
dummy_crise <- ts(dummy_crise,start = 1970,end = 2013,frequency = 1)
4 - Análise exploratória.
scatter.smooth(RGDP_FIJI,RGDP_AUS)
library(dygraphs)
dygraph(RGDP_FIJI)## Faz sentido fazermos um monthplot??
## monthplot(RGDP_FIJI)Vamos para a modelagem
PIB de Fiji e PIB da Austrália
Vamos começar tentando explicar o PIB de Fiji utilizando o PIB da Austrália, um parceiro comercial importante do país e que envia vários turistas para FIJI.
reg1<-lm(log(RGDP_FIJI) ~ log(RGDP_AUS))
summary(reg1)##
## Call:
## lm(formula = log(RGDP_FIJI) ~ log(RGDP_AUS))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12646 -0.02858 0.01698 0.03595 0.08670
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.19520 0.12383 -33.88 <2e-16 ***
## log(RGDP_AUS) 0.80096 0.01854 43.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05015 on 42 degrees of freedom
## Multiple R-squared: 0.978, Adjusted R-squared: 0.9775
## F-statistic: 1867 on 1 and 42 DF, p-value: < 2.2e-16## Teste do resíduo
plot(reg1$residuals,type="l")
abline(h=mean(reg1$residuals), col="red",lwd=1)
## Teste de autocorrelação
require(forecast)## Warning in sample.int(.Machine$integer.max - 1L, 1L): '.Random.seed' is not
## an integer vector but of type 'NULL', so ignoredcheckresiduals(reg1)
##
## Breusch-Godfrey test for serial correlation of order up to 10
##
## data: Residuals
## LM test = 17.909, df = 10, p-value = 0.05652require(lmtest)
dwtest(reg1)##
## Durbin-Watson test
##
## data: reg1
## DW = 0.68718, p-value = 9.965e-08
## alternative hypothesis: true autocorrelation is greater than 0Comentários:
1 - há problemas claros nos resíduos porque os resíduos apresentam um padrão de autocorrelação positiva;
2 - Essencialmente, a autocorrelação positiva está sugerindo o viés variável omitida;
3 - Além disso, temos um pequeno DW e R2 perto de 1 [regressão espúria].PIB de Fiji, PIB da Austrália e produção de açucar
reg2<-lm(log(RGDP_FIJI) ~ log(RGDP_AUS) + log(SUGAROUTPUT_FIJI))
summary(reg2)##
## Call:
## lm(formula = log(RGDP_FIJI) ~ log(RGDP_AUS) + log(SUGAROUTPUT_FIJI))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.038984 -0.013051 0.002852 0.012604 0.028664
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.92406 0.36643 -13.438 8.79e-12 ***
## log(RGDP_AUS) 0.78389 0.02889 27.132 < 2e-16 ***
## log(SUGAROUTPUT_FIJI) 0.10769 0.02305 4.671 0.000131 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01953 on 21 degrees of freedom
## (20 observations deleted due to missingness)
## Multiple R-squared: 0.9868, Adjusted R-squared: 0.9855
## F-statistic: 784.9 on 2 and 21 DF, p-value: < 2.2e-16## Teste do resíduo
plot(reg2$residuals,type="l")
abline(h=mean(reg2$residuals), col="red",lwd=1)
## Teste de autocorrelação
require(forecast)
checkresiduals(reg2)
##
## Breusch-Godfrey test for serial correlation of order up to 10
##
## data: Residuals
## LM test = 15.294, df = 10, p-value = 0.1217require(lmtest)
dwtest(reg2)##
## Durbin-Watson test
##
## data: reg2
## DW = 1.174, p-value = 0.006431
## alternative hypothesis: true autocorrelation is greater than 0Comentários:
1 - os resíduos estão começando a parecer um pouco mais aleatórios;
2 - No entanto, existem áreas problemáticas claras (20). E especificamente, eu vou apontar você para esta área aqui, o que na verdade é o período correspondente à Grande Recessão - ou a Grande Crise Financeira que ocorreu nos Estados Unidos;
3 - Outra possibilidade é o embargo europeu [vamos observar o gráfico abaixo].require(dygraphs)
dygraph(RGDP_FIJI)PIB de Fiji, PIB da Austrália, produção de açucar e uma dummy para a crise
O que faremos? Vamos criar uma variável dummy que nos permita modelar esse período atípico da nossa série de tempo.
reg3<-lm(log(RGDP_FIJI) ~ log(RGDP_AUS) + log(SUGAROUTPUT_FIJI) + dummy_crise)
summary(reg3)##
## Call:
## lm(formula = log(RGDP_FIJI) ~ log(RGDP_AUS) + log(SUGAROUTPUT_FIJI) +
## dummy_crise)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.02693 -0.01039 0.00277 0.01150 0.02932
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.96624 0.32784 -15.149 2.01e-12 ***
## log(RGDP_AUS) 0.79766 0.02639 30.225 < 2e-16 ***
## log(SUGAROUTPUT_FIJI) 0.10147 0.02075 4.891 8.83e-05 ***
## dummy_crise -0.02721 0.01084 -2.511 0.0207 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01745 on 20 degrees of freedom
## (20 observations deleted due to missingness)
## Multiple R-squared: 0.99, Adjusted R-squared: 0.9885
## F-statistic: 657.5 on 3 and 20 DF, p-value: < 2.2e-16## Teste do resíduo
plot(reg3$residuals,type="l")
abline(h=mean(reg3$residuals), col="red",lwd=1)
## Teste de autocorrelação
require(forecast)
checkresiduals(reg3)
##
## Breusch-Godfrey test for serial correlation of order up to 10
##
## data: Residuals
## LM test = 13.512, df = 10, p-value = 0.1964require(lmtest)
dwtest(reg3)##
## Durbin-Watson test
##
## data: reg3
## DW = 1.656, p-value = 0.08913
## alternative hypothesis: true autocorrelation is greater than 0Comentários:
1 - a variável dummy é significante e ajudou a modelar o período atípico ocorrido em nossa série temporal.PIB de Fiji: inserindo lag variables no modelo
Agora, como precisamos de variáveis de atraso, devemos instalar o pacote “dynlm”:
# install.packages("dynlm")
library(dynlm)
reg4 <- dynlm(log(RGDP_FIJI) ~ L(log(RGDP_FIJI),1) + log(RGDP_AUS) + L(log(RGDP_AUS),1) + log(SUGAROUTPUT_FIJI) + dummy_crise)
summary(reg4)##
## Time series regression with "ts" data:
## Start = 1990, End = 2013
##
## Call:
## dynlm(formula = log(RGDP_FIJI) ~ L(log(RGDP_FIJI), 1) + log(RGDP_AUS) +
## L(log(RGDP_AUS), 1) + log(SUGAROUTPUT_FIJI) + dummy_crise)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0237370 -0.0114449 -0.0000327 0.0108546 0.0241914
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.30599 0.65624 -6.562 3.64e-06 ***
## L(log(RGDP_FIJI), 1) 0.10749 0.13005 0.827 0.419330
## log(RGDP_AUS) 1.39504 0.28195 4.948 0.000104 ***
## L(log(RGDP_AUS), 1) -0.69496 0.26321 -2.640 0.016623 *
## log(SUGAROUTPUT_FIJI) 0.08271 0.02065 4.005 0.000830 ***
## dummy_crise -0.02311 0.00964 -2.398 0.027562 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01534 on 18 degrees of freedom
## (19 observations deleted due to missingness)
## Multiple R-squared: 0.993, Adjusted R-squared: 0.9911
## F-statistic: 512.4 on 5 and 18 DF, p-value: < 2.2e-16## Teste do resíduo
plot(reg4$residuals,type="l")
abline(h=mean(reg4$residuals), col="red",lwd=1)
## Teste de autocorrelação
require(forecast)
checkresiduals(reg4)
##
## Breusch-Godfrey test for serial correlation of order up to 10
##
## data: Residuals
## LM test = 15.142, df = 10, p-value = 0.127require(lmtest)
dwtest(reg4)##
## Durbin-Watson test
##
## data: reg4
## DW = 2.3243, p-value = 0.5943
## alternative hypothesis: true autocorrelation is greater than 0
Vamos ver a previsão dentro da amostra:
ts.plot(window(RGDP_FIJI,start=1990),exp(predict(reg4)),col=c(1,2))
PIB de Fiji: previsão fora da amostra
Comentários:
1 - Para fazer a previsão fora da amostra nós precisamos, prever as variáveis independentes [precisamos criar uma nova base de dados - x];
2 - Como nós não temos relações simultâneas entre as variáveis, não é necessário estimar a re4lação de longo prazo no modelo [A drunk and her dog](https://github.com/pedrocostaferreira/SeriesTemporais/blob/master/Aula9_ADL%2BECM_a-drunk-and-her-dog.pdf).
Observe que temos no nosso modelo como variáveis explicativas o PIB da Australia e a produçã de açucar de FIJI, dessa forma, precisamos criar uma base com essas variáveis para a previsão.
# dados para a previsão
x <- ts(Fiji_dataset[-nrow(Fiji_dataset), c("RGDP_AUS","SUGAROUTPUT_FIJI")], start = 1970, freq = 1)Agora, eu quero prever mais de uma série temporal simultaneamente. Para isso, estou usando as funções “lapply” e “auto.arima” para estimar automaticamente modelos ARIMA e fazer previsões para o PIB da Australia e a produção de açucar de FIJI.
library(forecast)
models <- lapply(x, FUN = forecast::auto.arima, allowdrift = F)
forecasts <- lapply(models, FUN = forecast::forecast, h = 4, level = 95)
RGDP_AUS_fcst <- forecasts$RGDP_AUS
SUGAROUTPUT_FIJI_fcst <- forecasts$SUGAROUTPUT_FIJIAgora, eu vou usar as informações obtidas acima para fazer a previsão 4 passos à frente:
% Previsão 1 passo à frente
# 2014 forecast
newdata <- matrix(c(1, # intercept
log(RGDP_FIJI[length(RGDP_FIJI)]), # log(RGDP_FIJI)_(t-1)
log(RGDP_AUS_fcst$mean[1]), # log(RGDP_AUS)
log(RGDP_AUS[length(RGDP_AUS)]), # log(RGDP_AUS)_(t-1)
log(SUGAROUTPUT_FIJI_fcst$mean[1]), # log(SUGAROUTPUT_FIJI)
0), # dummy_crise
nrow = 6, ncol = 1)
prev2014 <- exp(reg4$coefficients %*% newdata)% Previsão 2 passos à frente
# 2015 forecast
newdata <- matrix(c(1, # intercept
log(prev2014), # log(RGDP_FIJI)_(t-1)
log(RGDP_AUS_fcst$mean[2]), # log(RGDP_AUS)
log(RGDP_AUS_fcst$mean[1]), # log(RGDP_AUS)_(t-1)
log(SUGAROUTPUT_FIJI_fcst$mean[2]), # log(SUGAROUTPUT_FIJI)
0), # dummy_crise
nrow = 6, ncol = 1)
prev2015 <- exp(reg4$coefficients %*% newdata)% Previsão 3 passos à frente
# 2016 forecast
newdata <- matrix(c(1, # intercept
log(prev2015), # log(RGDP_FIJI)_(t-1)
log(RGDP_AUS_fcst$mean[3]), # log(RGDP_AUS)
log(RGDP_AUS_fcst$mean[2]), # log(RGDP_AUS)_(t-1)
log(SUGAROUTPUT_FIJI_fcst$mean[3]), # log(SUGAROUTPUT_FIJI)
0), # dummy_crise
nrow = 6, ncol = 1)
prev2016 <- exp(reg4$coefficients %*% newdata)% Previsão 4 passos à frente
# 2017 forecast
newdata <- matrix(c(1, # intercept
log(prev2016), # log(RGDP_FIJI)_(t-1)
log(RGDP_AUS_fcst$mean[4]), # log(RGDP_AUS)
log(RGDP_AUS_fcst$mean[3]), # log(RGDP_AUS)_(t-1)
log(SUGAROUTPUT_FIJI_fcst$mean[4]), # log(SUGAROUTPUT_FIJI)
0), # dummy_crise
nrow = 6, ncol = 1)
prev2017 <- exp(reg4$coefficients %*% newdata)% Todas as previsões fora da amostra
# all forecasts
all <- ts(c(exp(predict(reg4)), prev2014, prev2015, prev2016, prev2017), start = 1990, freq = 1)
ts.plot(window(RGDP_FIJI, start = 1990, freq = 1), all, col = 1:2)
all <- ts(c(prev2014, prev2015, prev2016, prev2017), start = 2014, freq = 1)
ts.plot(window(RGDP_FIJI, start = 1990, freq = 1), all, col = 1:2)
Prof. Dr. Pedro Costa Ferreira
Doutor em Engenharia Elétrica - (Decision Support Methods) e Mestre em Economia. Co-autor dos livros “Planejamento da Operação de Sistemas Hidrotérmicos no Brasil” e “Análise de Séries Temporais em R: um curso introdutório”. É o primeiro pesquisador da América Latina a ser recomendado pela empresa RStudio Inc. Atuou em projetos de Pesquisa e Desenvolvimento (P&D) no setor elétrico nas empresas Light S.A. (e.g. estudo de contingências judiciais), Cemig S.A, Duke Energy S.A, entre outras. Ministrou cursos de estatística e séries temporais na PUC-Rio e IBMEC e em empresas como o Operador Nacional do Setor Elétrico (ONS), Petrobras e CPFL S.A. Atualmente é professor de Econometria de Séries Temporais e Estatística, cientista chefe do Núcleo de Métodos Estatísticos e Computacionais (FGV|IBRE), coordenador dos cursos Economia Descomplicada (FGV|IDE) e Big Data e Data Science (FGV|IDE) e sócio-diretor da empresa Model Thinking Br (MTBr). É também revisor de importantes journals, como Energy Policy e Journal of Applied Statistics. Principais estudos são em modelos Econométricos, Incerteza Econômica, Preços, R software e Business Cycle.
contatos:
email: pedro.guilherme@fgv.br
website: pedrocostaferreira.github.io
GitHub: github.com/pedrocostaferreira
Linkedin: linkedin.com/pedro-costa-ferreira