1 Time Series

Conjunto de observações ordenadas no tempo:
* Uma possível trajetória de um processo aleatório
* Realização de um processo estocástico

Um TS pode ser contínua ou discreta. * Discreta: vendas diárias
* Contínua: Registro de venda de acordo com o instante que ocorre

Uma série é discreta quando se faz referencia ao tempo entre as observações e não à escala da variável. Tais séries contínuas podem ser discretizadas Barros et al. (2018). Via de regra TS são apresentadas como eventos discretos equi-espaçados.

A principal característica de séries temporais é a existência de alguma forma de dependência entre dados observados em tempos diferentes. O que limita a utilização de modelos estatísticos convencionais desenvolvidos para dados i.i.d.. Estes modelos podem ser utilizados para entender quais são os componentes envolvidos na série para assim tentar estimar valores futuros.

o domínio de TS pode ser paramétricos ou não paramétricos, entre modelos paramétricos podem ser citados modelos do tipo ARIMA, enquanto não paramétricos são considerados aqueles modelos de análise espectral a partir de decomposições de frequências.

Via de regra as TS Z(t) são modeladas sobre o parâmetro t fazendo referencia ao tempo, no entanto t pode ser dado em função de inúmeros parâmetros Morettin and Toloi (2020). Podendo as componentes serem r=3 [densidade, elasticidade, resistência] de dimenção p=3 [tempo, temperatura, altitude]; ou em um segundo exemplo casos de COVID-19 em duas cidades com r=1 [casos] e p=2 [tempo, cidade] Morettin and Toloi (2018). \[Z(t) = [Z_1(t),Z_2(t),Z_3(t)]\]

1.1 Processo Estocástico

Do ponto de vista teórico, uma TS é a realização de um processo estocástico. Dado um espaço de probabilidade \((\Omega, A, P)\), um processo estocástico seráa coleçãos de variáveis aleatórias no espaço amostral \(\Omega\), indexadas por um conjunto \(T\), denotado como:

\[\{x(t, \omega), t\in T\ e\ \omega \in \Omega\}\] Note que, entre todas as possibilidades, apenas observamos uma realização de tal processo estocástico.

Uma análise se TS visa:
Observar o mecanismo gerador da TS
Fazer previsões futuras de série (corto e longo prazo)
Descrever o comportamento da TS:
Graficos, diagramas de fispersão e histogramas, verificar tendências, ciclos e variabilidade sazonal.
Entender componentes constituintes:
Tendência: comportamento principal
Sazonalidade: padrão de repetiçõs por períodos diários, mensais, anuais, etc.

plot(AirPassengers)

ruído: diferença aleatória desacoplada de observações anteriores, da tandência e sazonalidade (média zero e comportamento estacionário - não correlacionado a x).

Algumas séries podem mudar sua tendenciaao longo do tempo, por exemplo deixando de ser linear para se tornar exponencial.

Simulação de um processo determinístico \(X_t = 1000 + 0.1t\)

n = 300
x = 1000+0.1*(1:n)
plot.ts(x)

Simulação de um processo de periodicidade 12 \(S_t = 20 sin((2\pi/12)t)\)

s = 20*sin(2*pi/12*(1:n))
plot.ts(s)

Somando termnos \(X_t + S_t\)

plot.ts(x+s)

Simulação de ruído branco $E_t N(0, 25) $

e = rnorm(n, 0.5)
plot.ts(e)

Somando termos \(X_t + S_t + E_t\)

z = x+s+e
plot.ts(z)

Decomposição

z = ts(z, start = 1990, frequency = 12)
plot(decompose(z))

Previsão

h        = 30
dec      = decompose(z)
trend    = dec$trend
sea      = dec$seasonal
e        = dec$random

tempo = 1:n
trend.lm = lm(trend~tempo)
trend.forec = trend.lm$coefficients[[1]] + trend.lm$coefficients[[2]]*((n+1):(n+h)) 
   c(trend, trend.forec) %>%ts(start=end(trend)+c(0,1), frequency = 12) %>%
   plot.ts()

SNAVE

library(forecast)
sazonal = snaive(sea, h= h)$mean

cbind(sea, sazonal) %>% plot.ts(plot.type = 's', col=c(1,2))

 previsao = trend.forec + sazonal
 cbind(z, previsao) %>% plot.ts(plot.type = 's', col=c(1,2))

Simulação de 100 pontos \(x_t = 5 + 0.35x_{t-1} + 0.6x_{t-6} + \epsilon_t\) com $ i.i.d.N(0,4) $

n = 100 
x = numeric(100)
x[1:6] = c(58.9,60.0,59.8,59.2,58.9,58.7)
for (i in 7:n) {
  x[i] = 5+0.35*x[i-1]+0.6*x[i-6]+rnorm(1,0,2)
}
x = ts(x, start = 2000, frequency = 6)
plot(x)

x.obs = window(x, end=time(x)[90])
x.teste =window(x, start=time(x)[91])
plot(x.obs)

plot(x.teste)

Bootstrap

n_sim = 200
x_matriz = matrix(NA, nrow = n, ncol = n_sim)

Observados

for (a in 1:n_sim) {
  x_matriz[1:90, a] = x.obs
}

Dados Simulados

for (b in 91:100) {
  x_matriz[b, ] = 5+0.35*x_matriz[b-1,]+0.6*x_matriz[b-6,]+rnorm(n_sim,0,2)
}
x_matriz = ts(x_matriz, start = 2000, frequency = 6)
plot(x_matriz, plot.type = 'single')

LI = LS = Pontual = numeric(10)
for (i in 1:10) {
  LI[i] = quantile(x_matriz[90+i, ], prob=0.025)
  LS[i] = quantile(x_matriz[90+i, ], prob=0.975)
  Pontual[i] = mean(x_matriz[90+i,])
}

prev_mat = cbind(LI, Pontual, LS)
colnames(prev_mat) = c("LI95", "Pontual", "LS95")
prev_mat = ts(prev_mat, start = time(x)[91], frequency = 6)
plot(x, ylim = c(55, 100))
points(prev_mat[, 1], col=3, type = 'l')
points(prev_mat[, 2], col=2, type = 'l')
points(prev_mat[, 3], col=3, type = 'l')

2 Média, Autocorrelação e Estacionariedade

2.1 Distribuição de Probabilidade em TS

Seja o processo ou trajetória \(\{X_t\}\) desdobrado como \(\{X_t, t=1,2,...,n\}\) uma sequência de v.a.s no espaço amostral \(\Omega\), pelo teoria de probabilidade e o teorema de Bayes sobre a fatoração da distribuição conjunta como o produto seguinte: \[P(X_1= x_1,...,X_n= x_n) = P(X_1= x_1)P(X_2= x_2|X_1= x_1),...,P(X_n= x_n|X_1= x_1,...,P(X_2= x_2|X_{n-1}= x_{n-1})\] Temos a função de probabilidade \(P(.)\) para processsos discretos e da mesma forma a f.d.p. para processos contínuos, onde:
Variável aleatória: \(x_t\)
Processo ou trajetória: \({x_t}\)
Distribuição de probabilidade: \(p(x_i)= P(X_i=x_i)\)
Distribição condicional: \(p_j(x_i)= P(X_i=x_i|X_i=x_i,...,X_j=x_j)\)
Distribuição conjunta: $ {i = 1}^{n} p{i-1}(x_i)$

2.2 Ruído Branco

  • Média Constante
  • Variância Constante
  • Sem correlacão temporal

Portanto não possui um padrão, sendo completamente aleatório

Simulando ruído branco com componentes autorregressivo, de diferenciação e de média movel (p, d, q) = (0, 0, 0).

whitenoise = arima.sim(model = list(order = c(0, 0, 0)), n = 50)
class(whitenoise)
## [1] "ts"
plot(whitenoise, col = "#333333", lwd = 2)

arima(whitenoise, order = c(0, 0, 0))
## 
## Call:
## arima(x = whitenoise, order = c(0, 0, 0))
## 
## Coefficients:
##       intercept
##         -0.0598
## s.e.     0.1247
## 
## sigma^2 estimated as 0.7781:  log likelihood = -64.67,  aic = 133.35

Um processo \({x_t}\) que é i.i.d. cpm distribuição normal \(N(0, \sigma^2)\) com média e variância constantes.

library(astsa)
## 
## Attaching package: 'astsa'
## The following object is masked from 'package:forecast':
## 
##     gas
plot(diff(globtemp), col = "blue")

2.3 Passeio aleatório

  • Sem variância ou média especificada

  • Dependência ao longo do tempo

  • Mudança ou incremento são ruídos brancos

Cada momento seque um processo aditino em que o instante adial é a adição de algum ruído ao valor do dado anterior.

\[T_{atual} = T_{atual-1} + e_{atual}\\ \Rightarrow T_1 = T_{0} + e_{1}, \,\, t.q.\,\, \overline{x} = 0\] Um segundo modelo de passeio aleatório chamado drift adiciona uma constante ao modelo

\[T_{atual} = c + T_{atual-1} + e_{atual} \\ \Rightarrow c + T_1 = T_{0} + e_{1}, \,\, t.q.\,\, \overline{x} = 0\]

passeioaleat = arima.sim(model = list(order = c(0, 1, 0)), n = 50)
ts.plot(passeioaleat, col = "green", lwd = 2)

arima(passeioaleat, order = c(0, 0, 0))
## 
## Call:
## arima(x = passeioaleat, order = c(0, 0, 0))
## 
## Coefficients:
##       intercept
##         -1.2988
## s.e.     0.2198
## 
## sigma^2 estimated as 2.463:  log likelihood = -95.36,  aic = 194.71
passeioaleat = diff(passeioaleat)
ts.plot(passeioaleat, col = "green", lwd = 2)

arima(passeioaleat, order = c(0, 0, 0))
## 
## Call:
## arima(x = passeioaleat, order = c(0, 0, 0))
## 
## Coefficients:
##       intercept
##         -0.0126
## s.e.     0.1497
## 
## sigma^2 estimated as 1.12:  log likelihood = -73.78,  aic = 151.56

Em um passeio aleatório inexistem componentes de tendência ou sazonalidade, cada termo é o resultado do termo anterior acrescido de um termo de erro (ruído branco). Nesse caso o erro não pussui média zero.
\[x_t = x_{t-1}+e_t\] Dessa forma quando se conhece \(x_0\):

\[x_1 = x_0+e_1 \to x_1 \sim N(x_0, \sigma^2) \\ x_2 = x_1+e_2 \to x_2|x_1 \sim N(x_1, \sigma^2) \\ .\\ .\\ x_n = x_{n-1}+e_n \to x_n|x_{n-1} \sim N(0, \sigma^2) \\\] Em que esta função possui o seguinte kernel: \[p(x_1,...,x_n) = \prod_{t=1}^{n} p_{t-1}(x_t) \\ \rightarrow \prod_{t=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}}exp\Bigg[-\frac{x_t-x_{t-1}}{2\sigma^2} \Bigg]\]

library(stats)
passeioaleat = arima.sim(model = list(order = c(0, 1, 0)), n = 50)
ts.plot(passeioaleat, col = "green", lwd = 2)

2.4 Correlacão

nn <- 100
tt1 <- sin(2 * pi * seq(nn)/12) - seq(nn)/50
tt2 <- cos(2 * pi * seq(nn)/12) - seq(nn)/10
plot(tt1, tt2, main = paste("Correlação", cor(tt1, tt2), sep = ": "))

t1 = rnorm(100, 5, 1)
t3 = t1*pi*runif(1, -25.0, 100)+runif(1, -25.0, 100)
plot(t1, t3, main = paste("Correlação", cor(t1, t3), sep = ": "))

plot(t1, tt1, main = paste("Correlação", cor(t1, tt1), sep = ": "))

2.5 Autocorrelação

A função de autocovariancia é definida por

\[\gamma(s, t) = cov(x_s, x_t)=E[(x_s-\mu_t)(x_t-\mu_t)]\] Variancia: \(\gamma(t, t)\)

A funcao de autocorrelacao como \[\rho(s, t) cor(x_s, x_t)=\frac{\gamma(s,t)}{\sqrt{\gamma(s, t)}\sqrt{\gamma(t, t)}}=\frac{\gamma(s,t)}{\sigma_s \sigma_t}\] Obs: \[\rho(s, t)=\rho(t, s) \\ -1 \leq \rho(s, t) \leq 1 \\ \rho(s, t)\ mede\ a\ relacao\ linear\ entre\ dois\ pontos\ da\ serie\]

tt <- seq(nn)
acf(tt)

plot.acf <- function(ACFobj) {
    rr <- ACFobj$acf[-1]
    kk <- length(rr)
    nn <- ACFobj$n.used
    plot(seq(kk), rr, type = "h", lwd = 2, yaxs = "i", xaxs = "i", 
        ylim = c(floor(min(rr)), 1), xlim = c(0, kk + 1), xlab = "Lag", 
        ylab = "Correlation", las = 1)
    abline(h = -1/nn + c(-2, 2)/sqrt(nn), lty = "dashed", col = "blue")
    abline(h = 0)
}

nn <- 100
tt <- seq(nn)
par(mfrow = c(1, 2))
plot.ts(tt, ylab = expression(italic(x[t])))
line.acf <- acf(tt, plot = FALSE)
plot.acf(line.acf)

tt <- sin(2 * pi * seq(nn)/12)
par(mfrow = c(1, 2))
plot.ts(tt, ylab = expression(italic(x[t])))
sine_acf <- acf(tt, plot = FALSE)
plot.acf(sine_acf)

tt <- sin(2 * pi * seq(nn)/12) - seq(nn)/100
par(mfrow = c(1, 2))
plot.ts(tt, ylab = expression(italic(x[t])))
sili_acf <- acf(tt, plot = FALSE)
plot.acf(sili_acf)

tt <- sin(2 * pi * seq(nn)/12) - seq(nn)/50
par(mfrow = c(1, 2))
plot.ts(tt, ylab = expression(italic(x[t])))
sili_acf <- acf(tt, plot = FALSE)
plot.acf(sili_acf)

tt <- sin(2 * pi * seq(nn)/12) - seq(nn)/10
par(mfrow = c(1, 2))
plot.ts(tt, ylab = expression(italic(x[t])))
sili_acf <- acf(tt, plot = FALSE)
plot.acf(sili_acf)

tt <- sin(2 * pi * seq(nn)/12) - seq(nn)/2
par(mfrow = c(1, 2))
plot.ts(tt, ylab = expression(italic(x[t])))
sili_acf <- acf(tt, plot = FALSE)
plot.acf(sili_acf)

2.6 Estacionariedade

As suposiciões iniciais de uma TS sao de que ela possui comportamento estacionário, perfazendo um trajeto ao longo do tempo com média constante; o que se traduz em um equilíbrio estável. Por mais estranho que pareça, a previsão de valores futuros requer tal equilíbrio, muito embora como senso comum se entanda que a série deva ser variável e não constante ao longo do tempo. Claro a maior parte das Ts apresentam uma não estacionariedade ou tendência - inclusive tendencia explosiva em casos exponenciais com o crescimento de uma colônia de bactérias ao longo do tempo Morettin and Toloi (2018).

Tais TS podem apresentar mudança de tendência podendo iniciar com uma certa inclinação moderada que pode sofrer alguma aceleração em um certo instante do tempo. Quando TS forem estacionárias e mesmo não estacionárias de maneira homogênea modelos ARIMA podem ser bastante úteis e de fácil aplicação, exceto naqueles casos exponenciais Morettin and Toloi (2018).

Como a implementação de modelos de projeção de TS requerem alguma estacionariedade, naqueles casos não estacionários são aplicados processos, como de diferenciação, para transformar as TS em modelos estacionários. Tal diferenciação consiste em obter diferenças sicessivas da TS original.

Estacionariedade Forte quando o processo for estritamente estacionário se a distribuição conjunta \(\{x_{t1},{x_{t2},...,{x_{tk} \}\) for a mesma com h defasagens \(\{x_{t1+h}, {x_{t2+h},...,{x_{tk+h}\}\).

Estacionariedade Fraca sera um processo com variancia finita com media constante \(\mu_t= E[x_t] = \mu\ \forall t=1,2,...\) e funcao de autocovariancias \(\gamma(s, t)\) dependendo apenas de s e t a partir da diferenca \(|s-t|\). Assim se o provesso for estacionario \(var[x_t]\) e \(E[x_t^ 2]\) serao constantes. Sendo s e t defasagens. Um processo sera estacionario de proimeira ordem se \(E[x_t]\) for constante e de segunda ordem se \(E[x_t]\) e \(E[x_t^2]\) forem constantes.

?Nile

Nile {datasets} R Documentation Flow of the River Nile Description Measurements of the annual flow of the river Nile at Aswan (formerly Assuan), 1871–1970, in 10^8 m^3, “with apparent changepoint near 1898” (Cobb(1978), Table 1, p.249).

Usage Nile Format A time series of length 100.

Source Durbin, J. and Koopman, S. J. (2001). Time Series Analysis by State Space Methods. Oxford University Press. http://www.ssfpack.com/DKbook.html

References Balke, N. S. (1993). Detecting level shifts in time series. Journal of Business and Economic Statistics, 11, 81–92. doi: 10.2307/1391308.

Cobb, G. W. (1978). The problem of the Nile: conditional solution to a change-point problem. Biometrika 65, 243–51. doi: 10.2307/2335202.

head(Nile)
## Time Series:
## Start = 1871 
## End = 1876 
## Frequency = 1 
## [1] 1120 1160  963 1210 1160 1160
RioNilo = Nile
tail(RioNilo)
## Time Series:
## Start = 1965 
## End = 1970 
## Frequency = 1 
## [1] 912 746 919 718 714 740
length(RioNilo)
## [1] 100
plot(RioNilo)

plot(RioNilo, type = "b", col = "green")

\(log()\) : Linearização
\(diff()\) : Remoção de tendência linear

2.7 Transformação

A aplicação de transformações busca obter séries com variancia controlada e efeito sazonal aditivo Morettin and Toloi (2018), para que a variação na amplitide da variância ao longo do processo não afete o modelo transformações como redução logarítimica podem ser aplicadas. Embora seja comum o uso de transformação sazonal, este não é ponto pacífico na literatura especializada. São comuns transformações como de Box-Cox, Bickel e Doksum, ou Hinkley.

Box-Cox: \[Z_t^{(\lambda)}= \begin{cases} \frac{z_t^{\lambda}-c}{\lambda},\ se\ \lambda \neq 0, \\ log(Z_t),\ se\ \lambda = 0, \end{cases} \]

Manly:

\[Z_t^{(\lambda)}= \begin{cases} \frac{e^{\lambda Z_t}-1}{\lambda},\ se\ \lambda \neq 0, \\ Z_t,\ se\ \lambda = 0, \end{cases} \] Bickel e Doksum:
\[Z_t^{(\lambda)}= \frac{|Z_t|^{\lambda}sgn(Y)-1}{\lambda}\]

2.8 ACF

So pode ser aplicado a TS estacionarias. A acima foi asumido que \(\gamma (s,t)=cov(x_s, x_t)\) dependendo apenas das defasagens s e t pela diferenca \(|s-t|\), isso afirma que a funcao de autocovariancias é constante entre diferentes pontos de mesma defasagem. Assim, dado um valor de defasagem h inteiro, logo \[\gamma(h) = \gamma(t+h, t)\] sera constante para todo 4=1,2,3,… Por isso a funcao de autocovariancias costuma ser chamada de funcao de defasagem \(\gamma(h)\).

Dado que \(x_1, x-2, ,...,x_n\) de uma TS constitui-se de um processo estacionario, logo $ {x} = _{t=1}^{n}$ sera o estimador da media \(\mu = \mu_t = E[x_t]\) constante. Assim como tambem a autocovariancia sera cosntante entre pontos de mesma defasagem h, a covariancioa amostral sera de ordem h por meio de: \(\hat{\gamma}(h)=\frac{1}{n}\sum_{t-i}^{n-h}(x_t-\bar{x})(x_{t+h}-\bar{x})\) como estimador de \(\gamma(h)=cov(x_t, x_{t-h})\).

Bem como a autocorrelacao amostral de ordem h \(\hat{\rho}(h)=\frac{\hat{\gamma}(h)}{\hat{\gamma(0)}}\) sera um estimador para \(\rho(h)=cpr(x_t, x_{t-h})\) onde \(\hat{\rho}(h)\) é tratado como função de autocolação (FAC) ou autocorrelation function (ACF).

2.9 CCF

So pode ser aplicado a TS estacionarias,

2.10 Diferenciação

A diferenciação é um dos modelos mais comumente aplicados. ASTSA: Applied Statistical Time Series analysis

library(astsa)
temperatura_global = globtemp
plot(temperatura_global, col = "red")

Diferenciação para obter uma série estacionária

plot(diff(temperatura_global))

par(mfrow = c(2,1))
plot(temperatura_global, col = "red")
plot(diff(temperatura_global), col = "blue")

3 Decomposição

A decomposição de uma TS pode ser de forma aditiva ou multiplicativa

\[X_t = T_t + S_t + E_t \\ X_t = T_t * S_t * E_t \\\] Tal decomposição visa antes de tudo decompor a TS em componentes mais simples de serem modelados. para obter uma TS sem tendência ou sem sazonalidade basta subrimir o termo correspondente.

Sem Tendência \[Z_t = X_t -T_t \\ X_t = X_t / T_t \\\] Sem Sazonalidade \[T_t = X_t - S_t \\Y_t = X_t / S_t \\\]

plot(decompose(AirPassengers, type = "a"))

plot(decompose(AirPassengers, type = "m"))

3.1 Remoção de Tendência

A decomposição parte da remoção de tendência.

3.1.1 Linear

Como visto pode possuir tendência linear ou polinomial, se TS aditiva ou multiplicativa. Dessa forma a tendência pode ser removida a removendo essa componente de linearidade ou polinomial.

ts = AirPassengers
tempo = time(ts)
fit = lm(ts~tempo)
plot(ts)
lines(x = tempo, y=fit$fitted.values, type = 'l', col="blue")

3.1.2 Polinomial

ts = Nile
time = as.numeric(time(ts))
out2 = lm(ts~poly(time, 2))
out6 = lm(ts~poly(time, 6))
plot(ts)
lines(x = time, y=out2$fitted.values, type = 'l', col="blue")
lines(x = time, y=out6$fitted.values, type = 'b', pch=18, col="red", lty=2)
legend(1880, 600, legend=c("Polinômio de dois Termos", "Polinômio de seis Termos"),
       col=c("blue", "red"), lty=1:2, cex=0.8)

3.2 Filtros para Remoção de Tendência

A maneira mais corriqueira de remoção de tendência é por meio do uso de filtros de suavização.

3.2.1 Médias Móveis

Dentre eles um bastante comum é o filtro de médias móveis de valores do passado filter com sides = 1, que como será visto é um componente de modelos ARMA.

\[SMA_T = \frac{1}{q} \sum_{i = 0}^{q} X_{t-i}\] Quando maior o número de períodos da média móvel, mas suave será a linha de tendência.

ts = AirPassengers
SMA6 = stats::filter(ts, rep(1/6, 6), sides = 1)
SMA12 = stats::filter(ts, rep(1/12, 12), sides = 1)
time = as.numeric(time(ts))
plot(ts)
lines(x = time, y=SMA6, type = 'l', col="blue")
lines(x = time, y=SMA12, type = 'b', pch=18, col="red", lty=2)
legend(1950, 500, legend=c("Média Móvel de 6 períodos", "Média Móvel de 12 períodos"),
       col=c("blue", "red"), lty=1:2, cex=0.8)

Também é possível realizar médias móveis simétricas em que são utilizados valores antes a após o ponto de medida, penalizando valores no inínio e ao fim da série.

Quando esta for simétrica e q impar: \[MM_T = \frac{1}{q}\sum_{i=-s}^{s} X_{t+1},\ onde\ s = (q-1)/2\]

Quando esta for simétrica e q par: \[MM_T = \frac{1}{q}\sum_{i=-s+1}^{s} X_{t+1},\ onde\ s = q/2\] Tomando uma média móvem simétrica de ordem 3: \[MM_T = \frac{1}{3}(X_{t-1}+X_{t}+X_{t+1})\] Tomando uma média móvem simétrica de ordem 4: Nesse caso por não ser simétrica ao redor de t o algoritmo utiliza um ponto a mais a frente (futuro). \[MM_T = \frac{1}{3}(X_{t-1}+X_{t}+X_{t+1}+X_{t+2})\]

MMS3 = stats::filter(ts, rep(1/3, 3), sides = 2)
MMS4 = stats::filter(ts, rep(1/4, 4), sides = 2)
plot(ts)
lines(x = time, y=MMS3, type = 'l', col="blue")
lines(x = time, y=MMS4, type = 'l', pch=18, col="red", lty=1)
legend(1950, 500, legend=c("Média Móvel de ordem 3", "Média Móvel de ordem 4"),
       col=c("blue", "red"), lty=1, cex=0.8)

3.2.1.1 Combinando MA de valores passados e simétricos

MMR5 = stats::filter(ts, rep(1/5, 5), sides = 1)
MMS5 = stats::filter(ts, rep(1/5, 5), sides = 2)
plot(ts)
lines(x = time, y=MMR5, type = 'l', col="blue")
lines(x = time, y=MMS5, type = 'l', pch=18, col="red", lty=1)
legend(1950, 500, legend=c("Média Móvel Retroativa de ordem 5", "Média Móvel Simétrica de ordem 5"),
       col=c("blue", "red"), lty=1, cex=0.8)

3.2.2 Filtro Linear

Nesse caso utilizando tewrmos (p, q) com coeficientes \(a_i\) como pesos que somados devem ser \(\sum a = 1\).

\[MM_l = \sum_{i=-p}^{q} a_iX_{t+i}\] Assim: \[MM_L = 0.2X_{t-i}+0.6X_{t}+0.2X_{t+i}\] Fornece:

MMl3 = stats::filter(ts, c(0.2, 0.6, 0.2), sides = 2)
MMl5 = stats::filter(ts, c(0.018, 0.2, 0.6, 0.2, 0.098), sides = 2)
plot(ts)
lines(x = time, y=MMl3, type = 'l', col="blue")
lines(x = time, y=MMl5, type = 'l', col="red")
legend(1950, 500, legend=c("Média Móvel de Filtro Linear de 3 termos", "Média Móvel de Filtro Linear de 5 termos"),
       col=c("blue", "red"), lty=1, cex=0.8)

3.3 Extração de Sazonalidade

3.3.1 Decomposição Aditiva

Supondo a decomposição de tendência aditiva por meio de \(TS_t = T_t + S_t + E+t\) a TS com tendência ajustada será dada por meio de:

\[Z_t = X_t-T_t = S_t+E_t\] Se a ordem de MA for 13 o ajuste é dados por uma MA simétrica de ordem 13 para considerar seis instantes antes a após o ponto de estimação:

ts = AirPassengers
autoplot(ts)

tendencia = stats::filter(ts, rep(1/13, 13), sides = 2)
z = ts-tendencia # TS sem tendência

plot(ts, main="TS Original", plot.type = "m", lwd=2)

plot(tendencia, main="Componente de Tendência", plot.type = "m", lwd=2)

plot(z, main="Modelo Ajustado sem Tendência", plot.type = "m", lwd=2)

Feito o ajuste de tendência (z) um novo processo de suavização é aplicado para extração da componente de sazonalidade. Nesse exeplo o processo consistirá em avaliar as médias de meses correspondentes onde S de janeiro corresponde à média dos Z dos mezes de janeiro de que se tem registro; o mesmo sendo aplicado aos demais meses.

ts = AirPassengers
ciclos = length(ts)/12
MediasSs = t(matrix(data = z, nrow = 12))
Ss = colMeans(MediasSs, na.rm = T)
Ss = ts(rep(Ss, ciclos), start=start(ts), frequency = frequency(ts))
plot(Ss, main="Componente de Sazonalidade", plot.type = "m", lwd=2)

erro = ts - tendencia - z
plot(cbind(ts, tendencia, z, erro), main="Modelo", plot.type = "m", lwd=2)

### Decomposição Aditiva
\[X_t = T_t*S_t*E_t\] Passa-se então ao ajuste da tendência da TS, obtendo a TS ajustada por \(Z_t = \frac{X_t}{T_t}\), obtendo dessa forma a série \(S_t\) de \(Z_t\) por suavização sazonal com resíduo \(E_t = \frac{X_t}{T_t*S_t}\).

5 Modelo MSTL

a função MSTL do pacote forecast automatiza a função STL buscando otimizar parametros para t.windown e s.windown, adota multiplas sazonalidades e aplica transformações Box-Cox.

ts = AirPassengers
library(forecast)
res = mstl(ts, lambda = NULL)
autoplot(res)

Box-Cox

#Box-Cox
res = mstl(ts, lambda = 'auto')
autoplot(res)

Série com dados intradiários de consumo de energia com dados fornecidos a cada meia hora.

autoplot(taylor)

res = mstl(taylor, lambda = NULL)
autoplot(res)

Note que a componente sazonal se decompõe em duas, uma TS sazonal com ciclo de 48 correspondendo aos 48 momentos do dia e uma segunda de 336 relativa aos ciclos que compõem a semana.

library(fpp2)
ts = calls
autoplot(ts)

res = mstl(ts, lambda = NULL)
autoplot(res)

res = mstl(ts, lambda = NULL, window(ts, end = 4))#primeiras 4 semanas apenas
autoplot(res)

Referências

Barros, Anna Carolina, Daiane Marcolino de Mattos, Ingrid Christyne de Oliveira, Pedro Guilherme Costa Ferreira, and Victor Eduardo Leite de Almeida Duca. 2018. Análise de Séries Temporais em R: Um curso introdutório. Edited by FGV IBRE. 1st ed. Rio de Janeiro: ELSEVIER.
Morettin, Pedro A, and Clélia M C Toloi. 2018. Análise de Séries Temporais: Modelos Lineares Univariados. 3rd ed. São Paulo: Blucher.
———. 2020. Análise de Séries Temporais: Modelos Multivariados e não lineares. Vol. 2. São Paulo: Blucher.