4  Regressão de séries temporais

Autor

Robson Bruno Dutra Pereira

4.1 Regressão linear simples

A regressão de séries temporais serve para prever uma série em função de outra ou mais relacionadas. Não é o método mais usual e sofisticado, uma vez que os métodos mais avançados consideram o padrão da própria série para modelagem e previsão. Entretanto, em algumas situações tal abordagem pode ser útil. Ademais, há possibilidade de combinar a regressão com ARIMA, conforme será abordado futuramente.

Seja um problema onde deseja-se prever uma série temporal contínua, \(y_1, y_2, ..., y_T\), em função de outra, \(x_1, x_2, ..., x_T\).

Na Figura 7.6 são plotadas as séries temporais anuais de área plantada em hectare e grãos produzidos em toneladas no Brasil de 1977 a 2022, Ipeadata. As séries apresentam tendência de crescimento não linear.

Figura 4.1: Produção de grãos no Brasil

Ao plotar um diagrama de dispersão entre as séries, pode-se observar uma correlação linear positiva alta entre estas, conforme Figura 4.2.

Figura 4.2: Área plantada versus Produção de grãos

Conforme observado na Figura 4.3, pode-se considerar em diversos casos a aproximação de uma função linear para tal relação.

Figura 4.3: Modelo de regressão linear simples para produção de grãos em função da área plantada

O modelo linear plotado em azul pode ser descrito conforme Equação abaixo, onde \(\beta_0\) é uma constante e \(\beta_1\) é um coeficiente linear.

\[ \hat{y}_t = \beta_0+\beta_{1}x_t \]

As observações da série dependente ou resposta podem ser descritas conforme segue, como sendo o valor estimado, \(\hat y_t\), adicionado de um termo de erro ou resíduo, \(\varepsilon_t\).

\[ \begin{aligned} y_t = \hat{y}_t + \varepsilon_t \\ y_t = \beta_0 + \beta_1x_t + \varepsilon_t \\ \end{aligned} \]

Considerando as \(T\) observações das séries, \((x_1, y_1), (x_2, y_2), ..., (x_T, y_T)\), pode-se pensar em um modelo que minimize os erros de previsão. Uma vez que o erro é normalmente distribuído, com média nula e variância \(\sigma_\varepsilon^2\), sendo os resíduos normalmente distribuídos, com média nula e variância igual a \(\sigma_\varepsilon^2\), \(\varepsilon \sim N(0,\sigma_\varepsilon^2)\), pode-se trabalhar a minimização da soma dos quadrados dos erros de previsão, \(\sum_{t=1}^{T}\varepsilon_t^2\). A Figura Figura 4.4 plota em linhas vermelhas verticais os resíduos. O modelo plotado minimiza a soma dos quadrados de tais erros.

Figura 4.4: Erros de previsão para o modelo de regressão linear para produção de grãos em função da área plantada

A análise à seguir expõe os coeficientes do modelo estimado com teste t para significância destes. Neste curso não será dada ênfase na inferência, mas na previsão. De forma simples um valor t com alta magnitude ou um p-valor (Pr(>|t|)) baixo indica risco baixo de rejeitar a hipótese nula de ausência de efeito da variável ou série independente (\(x_t\)).

Series: producao 
Model: TSLM 

Residuals:
   Min     1Q Median     3Q    Max 
-40729 -18733   1072  18188  35065 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.726e+05  1.248e+04  -13.84   <2e-16 ***
area_plantada  6.285e+00  2.561e-01   24.54   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21080 on 46 degrees of freedom
Multiple R-squared: 0.929,  Adjusted R-squared: 0.9275
F-statistic: 602.2 on 1 and 46 DF, p-value: < 2.22e-16

O modelo para prever a produção de grãos em função da área plantada fica:

\[ \hat{y}_t = -1,89 \times 10^{5} + 6,782x_t \]

4.2 Regressão linear múltipla

No caso de onde há múltiplas séries independentes ou regressoras de interesse, \(x_{t1}, x_{t2}, ..., x_{tk}\), onde \(k\) é o número de séries consideradas para estimar \(y_t\), pode-se considerar o modelo com um coeficiente linear associado a cada variável, isto é:

\[ \hat{y}_t = \beta_0 + \beta_1x_{t1} + \beta_2x_{t2} + \cdots + \beta_kx_{tk} = \beta_0 + \sum_{j=1}^{k}\beta_jx_{tj}, \]

ou de forma matricial com \(\mathbf{X}_{[T\times (k+1)]}\) e \(\mathbf{\beta}_{[(k+1) \times 1]}\):

\[ \begin{aligned} \hat{\mathbf{y}} = \mathbf{X}\mathbf{\beta} \end{aligned}, \]

onde a matrix \(\mathbf{X}\) contém uma coluna unitária para a constante e uma coluna para cada série independente:

\[ \mathbf{X} = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1k}\\ 1 & x_{21} & x_{22} & \cdots & x_{2k} \\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{T1} & x_{T2} & \cdots & x_{Tk} \\ \end{bmatrix}, e\\ \]

e \(\beta\) consite no vetor de coeficientes:

\[ \mathbf{\beta}^T = \begin{bmatrix} \beta_0 & \beta_1 & \cdots & \beta_k\\ \end{bmatrix}. \\ \]

Os valores observados da série \(\mathbf{y}\) podem ser recuperados somando os valores preditos e o erro.

\[ \begin{aligned} \mathbf{y} = \mathbf{X}\mathbf{\beta} + \mathbf{\varepsilon} \end{aligned} \]

Tomando tal notação, a soma dos quadrados dos erros pode ser descrita como \(\sum_{t=1}^{T}\varepsilon_t^2 = \mathbf{\varepsilon}^T\mathbf{\varepsilon}\). Desenvolvendo tal expressão tem-se:

\[ \begin{aligned} L(\mathbf{\beta}) = \mathbf{\varepsilon}^T\mathbf{\varepsilon} = (\mathbf{y} - \mathbf{X}\mathbf{\beta})^T(\mathbf{y} - \mathbf{X}\mathbf{\beta}) \\ \mathbf{y}^T\mathbf{y} - 2\mathbf{\beta}^T\mathbf{X}^T\mathbf{y} + \mathbf{\beta}^T\mathbf{X}^T\mathbf{X}\mathbf{\beta} \end{aligned} \]

Para minimizar \(L\) em relação à estimativa de \(\mathbf{\beta}\), pode-se diferenciar tal quantidade em relação à \(\mathbf{\beta}\) e igualar a zero:

\[ \begin{aligned} \frac{\partial L}{\partial \mathbf{\beta}} = -2\mathbf{X}^T\mathbf{y} + 2\mathbf{X}^T\mathbf{X}\mathbf{\beta} = 0 \\ \hat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}(\mathbf{X}^T\mathbf{y}) \end{aligned} \]

A estimativa dos coeficientes \(\hat\beta\) obtida minimiza a soma dos quadrados dos resíduos, \(\sum \varepsilon_i^2\), sendo portanto as estimativas de mínimos quadrados.

Seja a série multivariada de índices de produção de bens de capital, de consumo duráveis, intermediários e de consumo não duráveis, plotada na Figura 4.5.

Figura 4.5: Bens de consumo duráveis, intermediários e de consumo não duráveis

Observa-se na Figura 4.6 diagramas de dispersão e correlações aos pares, ficando clara a relação linear entre tais índices.

Figura 4.6: Bens de consumo duráveis, intermediários e de consumo não duráveis

Uma vez que os bens de consumo duráveis dependem dos bens de capital, os quais consistem nos ativos de empresas para produzir produtos ou serviços, e dos bens intermediários, supõe-se a possibilidade de obter um modelo de regressão múltipla desta série em função daquelas. A seguir apresenta-se o modelo de regressão múltipla estimado para a prever o índice de bens de consumos duráveis \(y_t\) em função dos de capital, \(x_{t1}\), e de produtos intermediários, \(x_{t2}\). Ambos apresentaram significância estatística.

Series: consumo_duraveis 
Model: TSLM 

Residuals:
    Min      1Q  Median      3Q     Max 
-38.361  -7.949   1.382   8.824  33.937 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -60.12405    9.55694  -6.291 1.28e-09 ***
capital          0.73046    0.05886  12.409  < 2e-16 ***
intermediarios   1.09143    0.12334   8.849  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.04 on 268 degrees of freedom
Multiple R-squared: 0.7734, Adjusted R-squared: 0.7717
F-statistic: 457.4 on 2 and 268 DF, p-value: < 2.22e-16

O modelo obtido pode ser escrito confirme segue:

\[ \hat y_t = -60,12 + 0,73x_{t1} + 1,09x_{t2} \]

4.3 Desempenho do modelo de regressão

Uma forma de medir o ajuste do modelo obtido aos dados seria a partir do cálculo do coeficiente de determinação múltipla, \(R^2\), conforme segue,

\[ \begin{align} R^2 &= 1- SS_{E}/SS_T \\ R^2 &= 1- \frac{\sum_{t=1}^{T}(y_t-\hat{y}_t)^2}{\sum_{i=t}^{T}(y_y-\overline{y})^2} = 1- \frac{\sum_{t=1}^{T}\varepsilon_t^2}{\sum_{i=t}^{T}(y_y-\overline{y})^2}, \end{align} \]

onde \(SS_E\) consiste na soma dos quadrados dos erros e \(SS_T\) consiste na soma dos quadrados total, ou no numerador da variância da série a ser predita.

O \(R^2\) sempre aumenta ao se adicionar novos termos no modelo, mesmo se estes não são significativos, uma vez que não leva em conta os graus de liberdade no cálculo. O coeficiente de determinação múltipla ajustado, \(R^2_{aj}\), é uma métrica mais interessante, já que considera as médias dos quadrados ao invés das somas dos quadrados dos erros, conforme segue:

\[ \begin{align} R^2_{aj} &= 1- MS_{E}/MS_T \\ R^2_{aj} &= 1- \frac{\sum_{t=1}^{T}(y_t-\hat{y}_t)^2/(T-K-1)}{\sum_{t=1}^{T}(y_t-\overline{y})^2/(T-1)}, \end{align} \]

onde \(K\) é o número de termos no modelo, \(MS_E\) consiste média dos quadrados dos erros, que consiste na estimativa de variância dos resíduos, enquanto \(MS_T\) consiste na média dos quadrados total, ou na variância da série a ser predita.

Para o modelo de regressão simples para produção de grãos em função da área plantada foi obtido um ajuste de mais de 91%, garantindo tal percentual de explicação da variabilidade da série de produção em função da área plantada. A Figura 4.7 plota a série de produção de grãos observada e a aproximada pelo modelo de regressão.

Figura 4.7: Produção de grãos

A Figura 4.8 plota a série do índice de produção de bens de consumo duráveis observada e a aproximada pelo modelo de regressão múltipla.

Figura 4.8: Bens de consumo duráveis

4.4 Diagnóstico dos resíduos

A Figura 4.9 apresenta alguns gráficos dos resíduos do modelo de produção de grãos em função da área plantada. Pode-se observar que o primeiro gráfico apresenta tendências em alguns instantes, com padrão não aleatório, indicando presença de autocorrelação. O correlograma confirma presença de autocorrelação nos resíduos até a defasagem de ordem 4. Logo, há informação importante não capturada pelo modelo de regressão.

Figura 4.9: Resíduos do modelo de produão de grãos em função da área plantada

4.5 Preditores úteis

Muitas séries temporais costumam apresentar tendência. Um modelo linear simples pode considerar o tempo, \(t\), como preditor, de forma a capturar a tendência.

\[ \hat{y}_t = \beta_0+\beta_{1}t \]

Pode-se também considerar variáveis dicotômicas ou binárias (dummy) para considerar a sazonalidade, feriados, ou algum evento especial. Por exemplo, para sazonalidade anual com uma série de frequência quadrimestral, pode-se considerar como variávis dummy as expostas na Tabela Tabela 4.1.

Tabela 4.1: Criação de variáveis dummy
Quadrimestre \(d_{1,t}\) \(d_{2,t}\) \(d_{3,t}\) \(d_{4,t}\)
Q1 1 0 0 0
Q2 0 1 0 0
Q3 0 0 1 0
Q4 0 0 0 1
Q1 1 0 0 0
\(\vdots\) \(\vdots\) \(\vdots\) \(\vdots\) \(\vdots\)

A Figura 4.10 exibe a série temporal de temperatura de São João del-Rei para os dias 15 a 28 de junho com frequência horária.

Figura 4.10: Temperatura em São João del-Rei nas duas últimas semanas de junho

A Figura 4.11 expõe o gráfico sazonal da mesma série. Observa-se sazonalidade diária e uma tendência de crescimento, especialmente da temperatura mínima, para os dias considerados.

Figura 4.11: Gráfico sazonal da temperatura em São João del-Rei

Um modelo de regressão considerando a tendência e sazonalidade para este caso deve considerar como variáveis dummy o dia, de forma a capturar a variação hora a hora. Logo, neste caso, 24 variáveis dummy são criadas. Observa-se que além da tendência a maior parte de tais variáveis apresenta significância e o modelo contempla mais de 95% de variabilidade da série.

Series: Temp 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-3.24262 -0.76670 -0.08595  0.83706  2.98696 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   14.8450939  0.3398667  43.679  < 2e-16 ***
trend()        0.0092571  0.0006713  13.790  < 2e-16 ***
season()day2  -0.9164000  0.4499590  -2.037 0.042534 *  
season()day3  -1.5470856  0.4499605  -3.438 0.000665 ***
season()day4  -2.2770295  0.4501793  -5.058 7.25e-07 ***
season()day5  -3.2077152  0.4501588  -7.126 7.25e-12 ***
season()day6  -3.9241152  0.4501393  -8.718  < 2e-16 ***
season()day7  -4.4405151  0.4501208  -9.865  < 2e-16 ***
season()day8  -5.0712008  0.4501032 -11.267  < 2e-16 ***
season()day9  -5.5018864  0.4500867 -12.224  < 2e-16 ***
season()day10 -6.0897150  0.4500712 -13.531  < 2e-16 ***
season()day11 -6.2989721  0.4500567 -13.996  < 2e-16 ***
season()day12 -5.1868006  0.4500432 -11.525  < 2e-16 ***
season()day13 -2.0460577  0.4500306  -4.546 7.82e-06 ***
season()day14  2.7518281  0.4500191   6.115 2.89e-09 ***
season()day15  6.7997138  0.4500086  15.110  < 2e-16 ***
season()day16  8.4261710  0.4499991  18.725  < 2e-16 ***
season()day17  9.2883425  0.4499906  20.641  < 2e-16 ***
season()day18  9.8719425  0.4499831  21.938  < 2e-16 ***
season()day19  9.9483997  0.4499766  22.109  < 2e-16 ***
season()day20  9.1105712  0.4499711  20.247  < 2e-16 ***
season()day21  7.9013141  0.4499666  17.560  < 2e-16 ***
season()day22  4.6491999  0.4499631  10.332  < 2e-16 ***
season()day23  3.0685142  0.4499605   6.820 4.76e-11 ***
season()day24  1.5949714  0.4499590   3.545 0.000453 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.19 on 311 degrees of freedom
Multiple R-squared: 0.9622, Adjusted R-squared: 0.9593
F-statistic: 330.3 on 24 and 311 DF, p-value: < 2.22e-16

Há uma tendência de aumento da temperatura nos dias considerados de 0,00888 °C por hora, ou 0,213°C por dia. Cada coeficiente das variáveis dummy expõe a diferença da respectiva hora em relação à primeira hora do dia. Por exemplo, às 6 da manhã a diferença de temperatura em relação a 1 da manhã é igual a -4.05°C. A Figura Figura 4.12 apresenta a série observada e a prevista pelo modelo de regressão.

Figura 4.12: Temperatura em São João del-Rei

A Figura 4.13 apresenta os valores previstos e observados plotados com as horas do dia separadas em cores distintas. Observa-se o excelente ajuste do modelo e a variação da temperatura segundo hora do dia.

Figura 4.13: Temperatura em São João del-Rei

4.6 Seleção de preditores

A Figura 4.14 exibe séries temporais de petróleo refinado e derivados produzidos no Brasil em \(m^3\), Dados estatísticos - Agência nacional de petróleo. O petróleo considera o volume nacional e importado.

Figura 4.14: Petróleo refinado e produção de derivados

A Figura 4.15 expõe um gráfico de correlação entre tais séries. As séries foram agrupadas segundo a magnitude e sinal das correlações. A série de petróleo refinado tem correlação positiva com asfalto, oléo diesel, outros não energéticos, querosene de avião, coque e gasolina e correlação negativa ou desprezível com as demais séries.

Figura 4.15: Correlação entre petróleo refinado e derivados

A seguir expõe-se um modelo de regressão múltipla do petróleo refinado em função do volume produzido dos derivados. Apenas GLP, coque, querosene de avião, querosene iluminante e solvente não apresentaram significância estatística.

Series: petroleo 
Model: TSLM 

Residuals:
    Min      1Q  Median      3Q     Max 
-606505  -83530    1027   87600  505807 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)             1.777e+06  2.207e+05   8.053 2.31e-14 ***
GLP                     3.063e-01  2.373e-01   1.291 0.197927    
asfalto                 1.161e+00  2.524e-01   4.600 6.39e-06 ***
coque                   1.437e+00  3.397e-01   4.229 3.18e-05 ***
gasolina                4.606e-01  8.168e-02   5.638 4.18e-08 ***
lubrificante            1.479e+00  9.627e-01   1.536 0.125627    
nafta                   7.169e-01  1.323e-01   5.419 1.29e-07 ***
oleo_combustivel        1.118e+00  6.843e-02  16.341  < 2e-16 ***
oleo_diesel             1.110e+00  4.892e-02  22.694  < 2e-16 ***
outros_nao_energeticos  9.948e-01  2.691e-01   3.697 0.000262 ***
parafina                8.141e-01  3.614e+00   0.225 0.821944    
querosene_aviao         1.458e+00  1.540e-01   9.472  < 2e-16 ***
querosene_iluminante   -1.496e+00  2.927e+00  -0.511 0.609642    
solvente                8.388e-01  8.403e-01   0.998 0.319019    
data                   -8.128e+01  9.855e+00  -8.248 6.28e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 158500 on 281 degrees of freedom
Multiple R-squared:  0.96,  Adjusted R-squared: 0.958
F-statistic: 481.7 on 14 and 281 DF, p-value: < 2.22e-16

Geralmente, em regressão múltipla é interessante selecionar um modelo com os melhores preditores, de forma a melhorar a capacidade preditiva do modelo. Além do \(R^2_{aj}\), outras métricas podem ser usadas, por exemplo o critério de informação de Akaike, AIC. O AIC considera o número de parâmetros e o erro do modelo, sendo esta uma métrica mais interessante para seleção de modelos.

\[ AIC = T\text{ log}\bigg(\frac{SS_E}{T}\bigg)+2(k+2) \]

A Tabela 4.2 apresenta os resultados de AIC, AIC corrigido e \(R^2_{aj}\) para o modelo completo e um modelo sem os parâmetros GLP, lubrificante, parafina, querosene iluminante e solvente. O modelo reduzido apresentou melhor resultado, uma vez que minimiza AIC e AICc.

Tabela 4.2: \(R^2_{aj}\), AIC e AIC corrigido para os modelos
model adj_r_squared AIC AICc
completo 0.9580020 7104.830 7106.780
reduzido 0.9576833 7102.288 7103.218

A forma correta de realizar a remoção de coeficientes, conforme pode-se presumir pelo exemplo anterior, não é considerando a significância. Ademais, não sugere-se a remoção manual de coeficientes. Recomenda-se o uso de um algoritmo para tal fim. O algoritmo stepwise com eliminação para trás segue os seguintes passos:

  1. Comece com o modelo completo.
  2. Remova um preditor por vez. Mantenha o modelo com melhor desempenho, considerando por exemplo o AIC.
  3. Repita o procedimento até encontrar o melhor modelo.

4.7 Previsão

A previsão à frente não é geralmente possível quando se considera outras séries como preditores uma vez que não se conhece os valores futuros das séries. Por exemplo, no último caso do petróleo e derivados, os derivados vêm depois, sendo mais útil para previsão do volume refinado um modelo que considere ou a decomposição ou outros métodos mais sofisticados ainda não abordados, os quais levam em conta a autocorrelação da série. Entretanto, nos casos de regressão onde se utiliza como preditores a tendência e a sazonalidade a partir de variáveis dummy, é possível realizar a previsão à frente. A Figura Figura 4.16 expõe a previsão de dois dias à frente para a temperatura de São joão del-Rei utilizando o modelo anteriormente obtido.

Figura 4.16: Previsão da temperatura em São João del-Rei

Apesar de não ser possível realizar previsões à frente para casos de regressão considerando outras séries temporais como preditores, é possível realizar nestes casos previsões baseadas em cenários. Seja o caso exposto anteriormente onde deseja-se prever o índice de produção de bens de consumo duráveis em função dos bens de capital e intermediários. A Figura 4.17 expõe a previsão para dois cenários, um de aumento e outro de decréscimo nos índices de produção de bens de capital e intermediários. Este tipo de análise pode ser útil para adiantar possíveis cenários e viabilizar ações de planejamento.

Figura 4.17: Previsão de cenários para produção de bens de capital

4.8 Regressão não linear

Transformações podem ser interessantes para tratar a não lineariedade na série a ser predita e nas preditoras. Há muitos casos que a tendência observada em uma série não é linear. Nestes casos podem ser considerada transformação polinomial na tendência.

Uma opção interessante em alguns casos é utilizar transformação logarítimica em ambas séries independentes e dependente ou em apenas uma destas. Para o caso simples, com transformação em ambas as séries, tem-se o seguinte modelo.

\[ \text{log} (y_t) = \beta_0+\beta_1\text{log} (x_t) \]

Pode-se pensar em considerar um termo quadrático ou de ordem maior para a tendência, de forma a contemplar uma tendência não linear. Entretanto não se recomenda tal estratégia, uma vez que na maioria dos casos a previsão resultante extrapola muito a realidade.

\[ \hat{y}_t = \beta_0+\beta_{1}t+\beta_{2}t^2 \]

Uma abordagem mais interessante é pensar em modelos lineares por partes, de forma a considerar modelos lineares distintos em cada parte do horizonte temporal da série, consistindo em um modelo não linear formado por peças lineares (piecewise regression). Por exemplo, um modelo com uma divisão no tempo, \(\tau\), ´pode ser descrito conforme segue:

\[ \hat{y}_t = \beta_0+\beta_{1}t+\beta_{2}(t-\tau)_+, \] onde:

\[ (t-\tau)_+ = \bigg\{ \begin{align} 0 \text{, se } t< \tau \\ t-\tau \text{, se } t\geq \tau \end{align} \]

A Figura 4.18 expõe a série temporal da evolução do recorde na prova de 100m nado peito na categoria masculina. Pode-se observar que de 1960 até 1980 o melhor tempo diminuiu com maior intensidade, passando para outro padrão de variação de 1980 a 2000 e, posteriormente, de 2000 até aqui com outra tendência.

Figura 4.18: Evolução do recorde na prova de nado peito, 100m, masculino

A Figura 4.19 expõe a mesma série com três modelos, um com transformação logarítimica na resposta, resultando em um modelo exponencial. Também foi considerado um modelo linear e outro linear por partes (piecewise). Foram consideradas duas partições, em 1977 e 2000. Pode-se observar que o terceiro caso parece interessante para ajustar-se melhor em cada momento considerado da série.

Figura 4.19: Evolução do record na prova de nado peito, 100m, masculino com modelos não lineares

4.9 Implementação em R

A seguir são apresentadas boa parte das implementações na linguagem R para obter os dados, gráficos e análises expostos no presente capítulo.

Carregando pacotes.

Séries temporais anuais de área plantada em hectare e grãos produzidos em toneladas.

# metadata("DEPAE_SAFRA")
# available_series()
graos_plan <- ipeadata("DEPAE_SAFRAAREA")
graos_prod <- ipeadata("DEPAE_SAFRA")
graos  <- tsibble(valor = c(graos_plan$value,
                            graos_prod$value),
                  data = yearmonth(c(graos_plan$date,
                                     graos_prod$date)),
                  id = c(rep("area_plantada", nrow(graos_plan)),
                         rep("producao", nrow(graos_prod))),
                  key=id, index=data)

id.labs <- c("área plantada [ha]", "produção [ton]")
names(id.labs) <- c("area_plantada", "producao")

Visualizando.

graos |> autoplot(valor) +
  # geom_line(col = "slategrey") +
  facet_wrap(nrow=2, ~ id, scales = "free_y",
             labeller = labeller(id = id.labs)) + 
  guides(colour = "none") +
  labs(y = "", x = "", title = "Produção de grãos no Brasil")

Diagrama de dispersão e correlação entre as séries.

df_graos <- graos |>
            pivot_wider(names_from=id,
                        values_from=valor)

ggpairs(df_graos, columns = 2:3)

Modelo de regressão linear simples para produção em função da área plantada.

fit_graos <- df_graos |>
  model(TSLM(producao ~ area_plantada))

report(fit_graos)

Série ajustada para produção de grãos a partir do modelo de regressão obtido.

augment(fit_graos) |>
  ggplot(aes(x = data)) +
  geom_line(aes(y = producao, colour = "Data")) +
  geom_line(aes(y = .fitted, colour = "Fitted")) +
  labs(y = "Grãos [ton]",
    title = "Produção de grãos"
  ) +
  scale_colour_manual(values=c(Data="grey",Fitted="slateblue")) +
  guides(colour = guide_legend(title = NULL))

Resíduos para o modelo de produção de grãos.

fit_graos |> gg_tsresiduals()

Série temporal de temperatura de São João del-Rei.

tempo_sjdr <- read.csv("INMET_SJDR_2024.csv",
                       header=T)

tempo_sjdr$Hora..UTC. <- tempo_sjdr$Hora..UTC./100

tempo_sjdr <- tempo_sjdr[,1:3]

tempo_sjdr <- tempo_sjdr |>
  mutate(Data = dmy(Data)) |>
  mutate(Data = make_datetime(year(Data), month(Data), day(Data), Hora..UTC.)) |>
  select(!Hora..UTC.) |>
  as_tsibble(index = Data) |>
  rename(Temp = Temp..Ins...C.)
data_inicial <- as.POSIXct("2024-06-15 00:00:00")

datas_especificas <- seq(from = data_inicial, by = "hour", length.out = 14*24)

tempo_sjdr_14_dias <- tempo_sjdr |>
  filter(Data %in% datas_especificas)

Visualizando.

tempo_sjdr_14_dias |> autoplot(Temp) +
  labs(x = "", y = "Temp [°C]", 
       title="Temperatura em São João del-Rei")

Gráfico sazonal.

tempo_sjdr_14_dias |> 
  gg_season(Temp, period = "day") +
  labs(x = "Hora", y = "Temp [°C]",
       title="Temperatura: Gráfico sazonal")
fit_tempo <- tempo_sjdr_14_dias |>
  model(TSLM(Temp ~ trend() + season()))
report(fit_tempo)

Série ajustada.

augment(fit_tempo) |>
  ggplot(aes(x = Data)) +
  geom_line(aes(y = Temp, colour = "Data")) +
  geom_line(aes(y = .fitted, colour = "Fitted")) +
  scale_colour_manual(
    values = c(Data = "black", Fitted = "orangered")
  ) +
  guides(colour = guide_legend(title = "Series")) +
  labs(x = "", y = "Temp [°C]",
       title="Temperatura em São João del-Rei")

Previsão do tempo.

fc_tempo <- forecast(fit_tempo)

fc_tempo |>
  autoplot(tempo_sjdr_14_dias) +
  labs(
    title = "Previsão do tempo para São João del-Rei",
    y = "Temp [°C]"
  )

Série de produção de petróleo e derivados.

petro <- read.csv("petroleo_e_derivados.csv", header=T)
petro <- petro |> select(-c(gasolina_aviao,
                            outros_energeticos))
petro_ts <- petro |>
  pivot_longer(cols = petroleo:solvente) |>
  mutate(data = as.Date(data, format="%m/%d/%Y"))

petro_ts <- petro_ts |> 
  as_tsibble(key = name, 
             index = data)
petro_ <- petro_ts |>
  pivot_wider(names_from=name,
              values_from=value)
petro_ts |> autoplot(value) +
  facet_wrap(nrow=14, ~ name, scales = "free_y") +
  guides(colour = "none")

Correlação entre as séries.

R <- cor(petro[,2:15])
corrplot(R, method = 'color', order = 'hclust', type = 'upper')

Modelo de regressão múltipla da série de produção de petróleo em função dos derivados.

fit_petro <- petro_ |> 
  model(TSLM(petroleo ~ .-petroleo-data))
  report(fit_petro)

Série temporal do record de 100m nado peito masculino.

peito <- read.csv("nado_peito_100m_M.csv", header=T)

peito <- peito[-c(1:7,15,18,21,34),]

peito_ts <- tsibble(recorde = peito$Time,
                    date = as.Date(peito$Date, 
                                    format = "%m/%d/%Y"),
                    index=date)
peito_ts |> autoplot(recorde) +
  labs(y="tempo [s]", x="", title="Evolução do record na prova de nado peito")

Modelos com tendência linear, exponencial e de regressão por partes.

knots <- as.Date(c("1977-08-18", "2000-06-15"), format="%Y-%m-%d")

fit_trends <- peito_ts |>
  model(
    linear = TSLM(recorde ~ trend()),
    exponencial = TSLM(log(recorde) ~ trend()),
    `linear por partes` = TSLM(recorde ~ trend(knots = knots))
  )
fc_trends <- fit_trends |> forecast(h = 365*3)

peito_ts |>
  autoplot(recorde) +
  geom_line(data = fitted(fit_trends),
            aes(y = .fitted, colour = .model)) +
  autolayer(fc_trends, alpha = 0.5, level = 95) +
  labs(y = "tempo [s]",
       title = "Recorde nado peito")

4.10 Exercícios propostos

  1. Faça um modelo de regressão para a série de produção de bens de consumo duráveis em função da série de bens de capital e da série de bens intermediários. Interprete o modelo obtido.

  2. Faça a análise dos resíduos para o modelo obtido. Interprete os resultados.

  3. Obtenha os valores ajustados e visualize-os junto com a série de bens de consumo duráveis.

  4. A seguir expõe-se o código para realizar a previsão baseada em cenários para o modelo de regressão para produção de grãos em função da área plantada. Foram considerados dois cenários: um de aumento na área plantada nos dois próximos anos com 82000 e 85000 hectares plantados; outro de queda na área plantada nos dois próximos anos com 75000 e 72000 hectares.

future_scenarios <- scenarios(
  Aumento = new_data(df_graos, 2) |>
    mutate(area_plantada = c(82000,85000)),
  Queda = new_data(df_graos, 2) |>
    mutate(area_plantada = c(75000, 72000)),
  names_to = "Cenários")
fc <- forecast(fit_graos, new_data = future_scenarios)

df_graos |>
  autoplot(producao) +
  autolayer(fc) +
  labs(title = "Previsão de cenários para produção de grãos",        y = "")

Faça uma previsão de cenários usando o modelo de regressão para a série de produção de bens de consumo duráveis. Considere três meses à frente e um cenários de queda e outro de aumento nos índices de produção de bens de capital e da série de bens intermediários. A seguir expõe-se como seria implementado um cenário para este caso de regressão múltipla.

Aumento = new_data(df_prod, 3) |>
    mutate(capital = c(100,110,120), 
           intermediarios = c(90,100,110))