API tools faq
paste
Advertisement
diogodg

Untitled

diogodg
Sep 25th, 2020
631
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 14.78 KB | None | 0 0
raw download clone embed print report
  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 ###
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!