diogodg

Untitled

Sep 25th, 2020
274
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. library(readxl)
  2. library(urca)
  3. library(tseries)
  4. library(forecast)
  5. library(lmtest)
  6. library(FinTS)
  7. library(ggplot2)
  8.  
  9. ## importação dos dados
  10. dados <- read_excel(
  11.   "C:/Users/55319/Documents/Materiais - Ciências Econômicas/2020_1/Econometria II/serieagriepeccom020.xls")
  12.  
  13. ## Dados em série temporal
  14. export <- ts(dados$Export, start = c(1996,1), end = c(2019,12), freq = 12)
  15.  
  16. ### ETAPA 1. IDENTIFICAÇÃO ###
  17. ## 1.1 Gráfico da série
  18. plot(export,
  19.        main = "Exportações do setor da agricultura e pecuária - US$ (milhões)\nPeriodicidade: mensal",
  20.        xlab = "Período", ylab = "", col = "blue", lwd = 2)
  21.  
  22. ## 1.2 Identicação da presença de sazonalidade
  23. # 1.2.1 Função monthplot: detectar sazonalidade
  24. monthplot(export, main = "Monthplot", ylab = "")
  25.  
  26. # 1.2.2 Função boxplot: detectar sazonalidade
  27. boxplot(export ~ cycle(export), main = "Boxplot", ylab = "", xlab = "")
  28.  
  29. # 1.2.3 Decomposição da serie temporal
  30. plot(stl(export, s.window="periodic"))
  31. plot(decompose(export)) +
  32.   labs(title = "Decomposição dos componentes da série temporal")
  33.  
  34. ## 1.3 FAC e FACP
  35. par(mfrow=c(1,1))
  36. acf(export,lag.max=36, main = "FAC", xlab = "Defasagem", ylab = "")
  37. pacf(export,lag.max=36, main = "FACP", xlab = "Defasagem", ylab = "")
  38.  
  39. ## 1.4 Gráfico da série em primeira diferença
  40. par(mfrow=c(1,1))
  41. ts.plot(diff(export, lag = 1, differences = 1), main = "Gráfico da série em primeira diferença",
  42.         xlab = "Período", ylab = "")
  43.  
  44. ## 1.4 Testes de raiz unitária ##
  45. ## 1.4.1 Teste ADF (H0: não estacionário)
  46. # Em nível
  47. summary(ur.df(export, type=c("trend"),lags=24, selectlags = "BIC"))
  48. summary(ur.df(export, type=c("drift"),lags=24, selectlags = "BIC"))
  49. summary(ur.df(export, type=c("none"),lags=24, selectlags = "BIC"))
  50. # Em primeira diferença
  51. summary(ur.df(diff(export), type=c("trend"),lags=24, selectlags = "BIC"))
  52. summary(ur.df(diff(export), type=c("drift"),lags=24, selectlags = "BIC"))
  53. summary(ur.df(diff(export), type=c("none"),lags=24, selectlags = "BIC"))
  54. # Em primeira diferença sazonal
  55. summary(ur.df(diff(export, lag = 12), type = "trend", lags = 24, selectlags = "BIC"))
  56. summary(ur.df(diff(export, lag = 12), type = "drift", lags = 24, selectlags = "BIC"))
  57. summary(ur.df(diff(export, lag = 12), type = "none", lags = 24, selectlags = "BIC"))
  58.  
  59. ## 1.4.2 Teste de PP (H0: não estacionário)
  60. # Em nível
  61. summary(ur.pp(export,type=c("Z-tau"), model=c("trend"), lags=c("short")))
  62. summary(ur.pp(export,type=c("Z-tau"), model=c("constant"), lags=c("short")))
  63. # Em primeira diferença
  64. summary(ur.pp(diff(export, lag = 12),type=c("Z-tau"), model=c("trend"), lags=c("short")))
  65. summary(ur.pp(diff(export, lag = 12),type=c("Z-tau"), model=c("constant"), lags=c("short")))
  66.  
  67. ## 1.5 Comparando a variável em nível e em primeira diferença
  68. par(mfrow=c(2,1))
  69. plot(export,
  70.      main = "Valor das Exportações",
  71.      xlab = "Período", ylab = "")
  72. plot(diff(export),
  73.      main = "Valor das Exportações\n(Em primeira diferença)",
  74.      xlab = "Período", ylab = "")
  75.  
  76. # 1.5.1 FAC e FACP em nível e primeira diferença
  77. par(mfrow=c(2,2))
  78. acf(export,lag.max=36, main = "FAC", xlab = "Defasagem", ylab = "")
  79. pacf(export,lag.max=36, main = "FACP", xlab = "Defasagem", ylab = "")
  80. acf(diff(export),lag.max=36, main = "FAC - Diff", xlab = "Defasagem", ylab = "")
  81. pacf(diff(export),lag.max=36, main = "FACP - Diff", xlab = "Defasagem", ylab = "")
  82.  
  83. ### ETAPA 2. ESTIMAÇÃO ###
  84. ## 2.1 Candidatos a melhor modelo SARIMA(p,d,q)x(P,D,Q)
  85. fit1=coeftest(Arima(export, order=c(1,0,0), seasonal=list(order=c(0,1,0), period=12))); fit1     ## Candidato a melhor modelo ##
  86. fit2=coeftest(Arima(export, order=c(1,0,0), seasonal=list(order=c(1,1,0), period=12))); fit2     ## Candidato a melhor modelo ##
  87. fit3=coeftest(Arima(export, order=c(1,0,0), seasonal=list(order=c(0,0,1), period=12))); fit3     ## Candidato a melhor modelo ##
  88.  
  89. fit4=coeftest(Arima(export, order=c(1,1,0), seasonal=list(order=c(0,1,0), period=12))); fit4     ## Candidato a melhor modelo ##
  90. fit5=coeftest(Arima(export, order=c(1,1,0), seasonal=list(order=c(1,1,0), period=12))); fit5     ## Candidato a melhor modelo ##
  91. fit6=coeftest(Arima(export, order=c(1,1,0), seasonal=list(order=c(0,1,1), period=12))); fit6     ## Candidato a melhor modelo ##
  92.  
  93. fit7=coeftest(Arima(export, order=c(0,1,1), seasonal=list(order=c(0,1,0), period=12))); fit7     ## Candidato a melhor modelo ##
  94. fit8=coeftest(Arima(export, order=c(0,1,1), seasonal=list(order=c(1,1,0), period=12))); fit8     ## Candidato a melhor modelo ##
  95. fit9=coeftest(Arima(export, order=c(0,1,1), seasonal=list(order=c(0,1,1), period=12))); fit9     ## Candidato a melhor modelo ##
  96.  
  97. fit10=coeftest(Arima(export, order=c(1,1,1), seasonal=list(order=c(0,1,0), period=12))); fit10     ## Candidato a melhor modelo ##
  98. fit11=coeftest(Arima(export, order=c(1,1,1), seasonal=list(order=c(1,1,0), period=12))); fit11    ## Candidato a melhor modelo ##
  99. fit12=coeftest(Arima(export, order=c(1,1,1), seasonal=list(order=c(0,1,1), period=12))); fit12    ## Candidato a melhor modelo ##
  100.  
  101. fit13=coeftest(Arima(export, order=c(2,0,0), seasonal=list(order=c(0,1,0), period=12))); fit13    ## Candidato a melhor modelo ##
  102. fit14=coeftest(Arima(export, order=c(2,0,0), seasonal=list(order=c(1,1,0), period=12))); fit14    ## Candidato a melhor modelo ##
  103.  
  104. fit15=coeftest(Arima(export, order=c(3,0,0), seasonal=list(order=c(0,1,1), period=12))); fit15    ## Candidato a melhor modelo* ##
  105. # Nota: * Fit 15 não possui todos coeficientes significativos, porém possui mais componentes AR - vamos testar os resíduos com esse modelo
  106.  
  107. auto.arima(export)
  108. fit16=coeftest(auto.arima(export)); fit16      ## Candidato a melhor modelo ##
  109.  
  110. ## 2.2 Critérios de informação
  111. mod1=Arima(export, order=c(1,0,0), seasonal=list(order=c(0,1,0), period=12)); mod1     ## Candidato a melhor modelo ##
  112. mod2=Arima(export, order=c(1,0,0), seasonal=list(order=c(1,1,0), period=12)); mod2     ## Candidato a melhor modelo ##
  113. mod3=Arima(export, order=c(1,0,0), seasonal=list(order=c(0,0,1), period=12)); mod3     ## Candidato a melhor modelo ##
  114.  
  115. mod4=Arima(export, order=c(1,1,0), seasonal=list(order=c(0,1,0), period=12)); mod4     ## Candidato a melhor modelo ##
  116. mod5=Arima(export, order=c(1,1,0), seasonal=list(order=c(1,1,0), period=12)); mod5     ## Candidato a melhor modelo ##
  117. mod6=Arima(export, order=c(1,1,0), seasonal=list(order=c(0,1,1), period=12)); mod6     ## Candidato a melhor modelo ##
  118.  
  119. mod7=Arima(export, order=c(0,1,1), seasonal=list(order=c(0,1,0), period=12)); mod7     ## Candidato a melhor modelo ##
  120. mod8=Arima(export, order=c(0,1,1), seasonal=list(order=c(1,1,0), period=12)); mod8     ## Candidato a melhor modelo ##
  121. mod9=Arima(export, order=c(0,1,1), seasonal=list(order=c(0,1,1), period=12)); mod9     ## Candidato a melhor modelo ##
  122.  
  123. mod10=Arima(export, order=c(1,1,1), seasonal=list(order=c(0,1,0), period=12)); mod10     ## Candidato a melhor modelo ##
  124. mod11=Arima(export, order=c(1,1,1), seasonal=list(order=c(1,1,0), period=12)); mod11    ## Candidato a melhor modelo ##
  125. mod12=Arima(export, order=c(1,1,1), seasonal=list(order=c(0,1,1), period=12)); mod12    ## Candidato a melhor modelo ##
  126.  
  127. mod13=Arima(export, order=c(2,0,0), seasonal=list(order=c(0,1,0), period=12)); mod13    ## Candidato a melhor modelo ##
  128. mod14=Arima(export, order=c(2,0,0), seasonal=list(order=c(1,1,0), period=12)); mod14    ## Candidato a melhor modelo ##
  129.  
  130. mod15=Arima(export, order=c(3,0,0), seasonal=list(order=c(0,1,1), period=12)); mod15    ## Candidato a melhor modelo ##
  131.  
  132. mod16=auto.arima(export); mod16
  133.  
  134. AIC(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,
  135.     mod9, mod10, mod11, mod12, mod13, mod14, mod15, mod16)
  136. BIC(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,
  137.     mod9, mod10, mod11, mod12, mod13, mod14, mod15, mod16)
  138. # Nota:3 melhores modelos selecionados + modelo 15
  139. # Nota: Vamos considerar critério BIC, pois tende a ser mais parcimonioso
  140.  
  141. ### 3.VERIFICAÇÃO ###
  142. # Nota: Vamos avaliar os modelos 16, 12, 9 e 15
  143.  
  144. ## 3.1 Teste de estabilidade - saber se as raizes do polinômio característico do modelo satisfazem a propriedade de estacionariedade
  145. autoplot(mod16)
  146. autoplot(mod12)
  147. autoplot(mod9)
  148. autoplot(mod15)
  149. # Nota: todos os modelos são estáveis
  150.  
  151. ## 3.2 Testes dos resíduos
  152. # 3.2.1 Diagnóstico da série temporal
  153. tsdiag(mod16)
  154. tsdiag(mod12)
  155. tsdiag(mod9)
  156. tsdiag(mod15)
  157.  
  158. # 3.2.2 Autocorrelação                   Teste de Ljung-Box| H0: os residuos sao iid (não há autocorrelação)
  159. res12<-residuals(mod12)
  160. Box.test(res12,lag=12,type="Ljung-Box")
  161. Box.test(res12,lag=24,type="Ljung-Box")
  162. Box.test(res12,lag=36,type="Ljung-Box")
  163. # Nota: mod 12 tem autocorrelação
  164.  
  165. res16<-residuals(mod16)
  166. Box.test(res16,lag=12,type="Ljung-Box")
  167. Box.test(res16,lag=24,type="Ljung-Box")
  168. Box.test(res16,lag=36,type="Ljung-Box")
  169. # Nota: mod16 tem autocorrelação
  170.  
  171. tsdiag(mod9)
  172. res9<-residuals(mod9)
  173. Box.test(res9,lag=12,type="Ljung-Box")
  174. Box.test(res9,lag=24,type="Ljung-Box")
  175. Box.test(res9,lag=36,type="Ljung-Box")
  176. # Nota: mod 9 tem autocorrelação
  177.  
  178. tsdiag(mod15)
  179. res15<-residuals(mod15)
  180. Box.test(res15,lag=12,type="Ljung-Box")
  181. Box.test(res15,lag=24,type="Ljung-Box")
  182. Box.test(res15,lag=36,type="Ljung-Box")
  183. # Nota: tem autocorrelação
  184.  
  185. # 3.2.2 Normalidade                      Teste de Jarque-Bera| H0: normalidade dos residuos (normalmente distribuídos)
  186. # Análise visual mod16
  187. par(mfrow=c(2,2))
  188. hist(res16, freq=F, ylab='Densidade', xlab='Resíduos', main='Resíduos')
  189. plot(density(res16, kernel = c("gaussian")), main="Resíduos")   #Função de densidade estimada
  190. qqnorm(res16, ylab='Quantis amostrais', xlab='Quantis teóricos', main='Quantil-Quantil')
  191. qqline(res16, col = "red")
  192. # Nota: próximo do centro se aproxima aproxima a distribuição normal, mas nas extremidades se afasta bastante
  193. # Teste Shapiro e Jarque Bera
  194. shapiro.test(res16)
  195. jarque.bera.test(res16)
  196. # Nota: ambos rejeitam a hipótese de normalidade dos erros
  197.  
  198. # Análise visual mod12
  199. par(mfrow=c(2,2))
  200. hist(res12, freq=F, ylab='Densidade', xlab='Resíduos', main='Resíduos')
  201. plot(density(res12, kernel = c("gaussian")), main="Resíduos")   #Função de densidade estimada
  202. qqnorm(res12, ylab='Quantis amostrais', xlab='Quantis teóricos', main='Quantil-Quantil')
  203. qqline(res12, col = "red")
  204. # Teste Shapiro e Jarque Bera
  205. shapiro.test(res12)
  206. jarque.bera.test(res12)
  207. # Nota: ambos rejeitam a hipótese de normalidade dos erros
  208.  
  209. # Análise visual mod9
  210. par(mfrow=c(2,2))
  211. hist(res9, freq=F, ylab='Densidade', xlab='Resíduos', main='Resíduos')
  212. plot(density(res9, kernel = c("gaussian")), main="Resíduos")   #Função de densidade estimada
  213. qqnorm(res9, ylab='Quantis amostrais', xlab='Quantis teóricos', main='Quantil-Quantil')
  214. qqline(res9, col = "red")
  215. # Teste Shapiro e Jarque Bera
  216. shapiro.test(res9)
  217. jarque.bera.test(res9)
  218. # Nota: ambos rejeitam a hipótese de normalidade dos erros
  219.  
  220. # Análise visual mod15
  221. par(mfrow=c(2,2))
  222. hist(res15, freq=F, ylab='Densidade', xlab='Resíduos', main='Resíduos')
  223. plot(density(res15, kernel = c("gaussian")), main="Resíduos")   #Função de densidade estimada
  224. qqnorm(res15, ylab='Quantis amostrais', xlab='Quantis teóricos', main='Quantil-Quantil')
  225. qqline(res15, col = "red")
  226. # Teste Shapiro e Jarque Bera
  227. shapiro.test(res15)
  228. jarque.bera.test(res15)
  229. # Nota: ambos rejeitam a hipótese de normalidade dos erros
  230.  
  231. # 3.2.3 Teste de heteroscedasticidade     Teste ARCH| H0: os residuos nao possuem efeitos auto-regressivos de heteroscedasticidade condicional
  232.  
  233. ArchTest(mod16$residuals,lags = 36)
  234. ArchTest(mod12$residuals,lags = 36)
  235. ArchTest(mod9$residuals,lags = 36)
  236. ArchTest(mod15$residuals,lags = 36)
  237.  
  238. ### ETAPA 4. PREVISÃO ###
  239. ## 4.1 Testes de acurária
  240. # 4.1.1 Teste de acurácia: dentro da amostra (amostra inteira como treino - desvios do modelo ajustaoo em relação aos dados observados - média dos desvios ao quadrado)
  241. # Nota: vamos comparar o RMSE (root mean square error)
  242. accuracy(mod16)
  243. accuracy(mod12)
  244. accuracy(mod9)
  245. accuracy(mod15)
  246. ## modelo 12 e 16 erraram menos
  247.  
  248. # 4.1.2 Teste de acurácia: fora da amostra (usando o ano de 2019 como teste - cortar um pedaço da série, ajustar na parte restante e gerar previsões para a parte descartada)
  249. export.test <- tail(export,12)                                   # definindo a série "teste"
  250. export.test
  251.  
  252. export.train <- head(export, length(export)-12)                    # definindo a série "treino"
  253. export.train
  254.  
  255. # Modelo 16
  256. mod16.train=auto.arima(export.train); mod16.train
  257.  
  258. fc_mod16.train <- forecast(mod16.train, h = 12) # projeção para 12 períodos para frente
  259. fc_mod16.train
  260. par(mfrow=c(1,1))
  261.  
  262. accuracy(fc_mod16.train$mean,export.test)
  263.  
  264. df_test16 <- ts(data.frame(cbind(fcst=fc_mod16.train$mean,obs=export.test)))
  265. df_test16
  266. autoplot(df_test16[,2], series = "Observado") +
  267.   autolayer(df_test16[,1], series = "Previsão") +
  268.   labs(title = "Valor das Exportações (Período: Jan/2019 - Dez/2019)",
  269.        x = "Períodos",
  270.        y = "",
  271.        color = "Previsão")
  272. # Nota: modelo 16 conseguiu captar razoalvelmente bem as variações e a tendência dos valores observados
  273.  
  274. # Modelo 12
  275. mod12.train=Arima(export.train, order=c(1,1,1), seasonal=list(order=c(0,1,1), period=12)); mod12.train
  276.  
  277. fc_mod12.train <- forecast(mod12.train, h = 12)
  278. fc_mod12.train
  279.  
  280. accuracy(fc_mod12.train$mean,export.test)
  281.  
  282. df_test12 <- ts(data.frame(cbind(fcst=fc_mod12.train$mean,obs=export.test)))
  283. autoplot(df_test12[,2], series = "Observado") +
  284.   autolayer(df_test12[,1], series = "Previsão") +
  285.   labs(title = "Valor das exportações (Período: Jan/2019 - Dez/2019)",
  286.        x = "Períodos",
  287.        y = "",
  288.        color = "Previsão")
  289.  
  290. # Modelo 15
  291. mod15.train=Arima(export.train, order=c(3,0,0), seasonal=list(order=c(0,1,1), period=12)); mod15.train
  292.  
  293. fc_mod15.train <- forecast(mod15.train, h = 12)
  294. fc_mod15.train
  295.  
  296. accuracy(fc_mod15.train$mean,export.test)
  297.  
  298. df_test15 <- ts(data.frame(cbind(fcst=fc_mod15.train$mean,obs=export.test)))
  299. autoplot(df_test15[,2], series = "Observado") +
  300.   autolayer(df_test15[,1], series = "Previsão") +
  301.   labs(title = "Valor das exportações (Período: Jan/2019 - Dez/2019)",
  302.        x = "Períodos",
  303.        y = "",
  304.        color = "Previsão")
  305.  
  306. # Comparando acurácias
  307. accuracy(fc_mod16.train$mean,export.test)
  308. accuracy(fc_mod12.train$mean,export.test)
  309. accuracy(fc_mod15.train$mean,export.test)
  310. # Nota: pelos critérios de acurácia, escolheremos o modelo 16
  311.  
  312. # 4.2 Gráficos
  313. # Modelo 16
  314. # 4.2.1 Gráfico "manual"
  315. fc_mod16 <- forecast(mod16, h = 12)
  316. se16 <- sqrt(mod16$sigma2)
  317.  
  318. par(mfrow=c(1,1))
  319. ts.plot(window(export, start=c(2018,1)),
  320.         fc_mod16$mean,
  321.         fc_mod16$mean+1.96*se16,
  322.         fc_mod16$mean-1.96*se16,col=c(1,2,2,2), lty=c(1,1,2,2))
  323.  
  324. # 4.2.2 Gráfico da função 'forecast'
  325. plot(fc_mod16)
  326.  
  327. plot(fc_mod16,
  328.      main = "Previsão das Exportações - MG para 2021",
  329.      xlab = "Período",
  330.      include = 24,
  331.      showgap = F,
  332.      fcol = "darkgreen",
  333.      flty = "dashed")
  334.  
  335. ### FIM ###
RAW Paste Data