Modelos de previsão de tendências e sazonalidade com séries temporais

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon R-posts.com, e gentilmente contribuiu para os R-blogueiros]. (Você pode relatar um problema sobre o conteúdo desta página aqui)


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

Os preços da gasolina sempre são um problema na Turquia; porque o povo turco gosta de dirigir para onde quer que vá, mas de qualquer maneira reclamam dos preços. Eu queria começar a cavar nos preços dos últimos sete anos e como eles foram. Eu usei gasolina sem chumbo preços de 95 octanas da Petrol Ofisi, que é uma empresa de distribuição de combustível na Turquia. Organizei os preços mensalmente entre 2013 e 2020 como a moeda turca, Lira (TL). O conjunto de dados está aqui

head(gasoline_df)
#        date gasoline
#1 2013-01-01     4.67
#2 2013-02-01     4.85
#3 2013-03-01     4.75
#4 2013-04-01     4.61
#5 2013-05-01     4.64
#6 2013-06-01     4.72

Antes de tudo, eu queria ver como os preços da gasolina foram durante o período e analisar as linhas de tendência ajustadas.

#Trend lines 
library(tidyverse)

ggplot(gasoline_df, aes(x = date, y = gasoline)) + geom_point() +
  stat_smooth(method = 'lm', aes(colour = 'linear'), se = FALSE) +
  stat_smooth(method = 'lm', formula = y ~ poly(x,2), aes(colour = 'quadratic'), se= FALSE) +
  stat_smooth(method = 'lm', formula = y ~ poly(x,3), aes(colour = 'cubic'), se = FALSE)+
  stat_smooth(data=exp.model.df, method = 'loess',aes(x,y,colour = 'exponential'), se = FALSE) 

Como a variável de resposta é medida por logarítmica no modelo exponencial, precisamos obter o expoente e criar um quadro de dados diferente para a linha de tendência exponencial.

#exponential trend model
exp.model <- lm(log(gasoline)~date,data = gasoline_df) 
exp.model.df <- data.frame(x=gasoline_df$date,
                           y=exp(fitted(exp.model)))

Como vimos acima, os modelos cúbico e exponencial quase se sobrepõem e parecem se ajustar melhor aos dados. No entanto, analisaremos cada modelo em detalhes.

Os modelos de tendência de previsão

A tendência linear; y_ {t}, o valor da série em um determinado momento, t, é descrito como:

y_ {t} =  beta_ {0} +  beta_ {1} t +  epsilon_ {t}

 beta_ {0} e  beta_ {1} são os coeficientes.

model_linear <- lm(data = gasoline_df,gasoline~date)

Acima, criamos uma variável de modelo para o modelo de tendência linear. Para comparar os modelos, precisamos extrair os coeficientes de determinação ajustados, que são usados ​​para comparar modelos de regressão com um número diferente de variáveis ​​explicativas, de cada modelo de tendência.

Ajustado R ^ 2 = 1- (1-R ^ 2) ( frac {n-1} {n-k-1})
n: tamanho da amostra
k: número de variáveis
R ^ 2: coeficiente de determinação, que é o quadrado da correlação dos valores de observação com os valores previstos: r ^ 2_ {y  hat {y}}

adj_r_squared_linear <- summary(model_linear) %>% 
  .$adj.r.squared %>% 
  round(4)

A tendência exponencial; diferentemente da tendência linear, permite que a série aumente a uma taxa crescente a cada período, é descrito como:
 ln (y_ {t}) =  beta_ {0} +  beta_ {1} t +  epsilon_ {t}.  ln (y_ {t}) é um logaritmo natural da variável de resposta.

Para fazer previsões no modelo ajustado, usamos a função exponencial como  hat {y_ {t}} = exp (b_ {0} + b_ {1} t + s_ {e} ^ 2/2) porque a variável dependente foi transformada por uma função logarítmica natural. Para  hat {y_ {t}} não estar abaixo do valor esperado de y {t}, adicionamos metade do quadrado do erro padrão residual. A fim de encontrar s_ {e}, executamos a função de resumo como abaixo. Além disso, podemos ver o nível de significância dos coeficientes e do modelo. Como podemos ver abaixo, como o valor p dos coeficientes e do modelo é menor que 0,05, eles são significativos no nível de significância de% 5.

summary(model_exponential)

#Call:
#lm(formula = log(gasoline) ~ date, data = gasoline_df)

#Residuals:
#      Min        1Q    Median        3Q       Max 
#-0.216180 -0.083077  0.009544  0.087953  0.151469 

#Coefficients:
#              Estimate Std. Error t value Pr(<|t|)    
#(Intercept) -1.0758581  0.2495848  -4.311  4.5e-05 ***
#date         0.0001606  0.0000147  10.928  < 2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Residual standard error: 0.09939 on 82 degrees of freedom
#Multiple R-squared:  0.5929,  Adjusted R-squared:  0.5879 
#F-statistic: 119.4 on 1 and 82 DF,  p-value: < 2.2e-16

Criamos uma variável de modelo exponencial.

model_exponential <- lm(data = gasoline_df,log(gasoline)~date)

Para se ajustar R ^ 2, primeiro transformamos a variável prevista.

library(purrr)

y_predicted <- model_exponential$fitted.values %>% 
  map_dbl(~exp(.+0.09939^2/2))

E então, executamos a fórmula mostrada acima para nos ajustar R ^ 2

r_squared_exponential <- (cor(y_predicted,gasoline_df$gasoline))^2

adj_r_squared_exponential <- 1-(1-r_squared_exponential)*
  ((nrow(gasoline_df)-1)/(nrow(gasoline_df)-1-1))

Às vezes, uma série temporal muda a direção em relação a vários motivos. Esse tipo de variável de destino é ajustado para modelos de tendência polinomial. Se houver uma mudança de direção, usamos o modelo de tendência quadrática.

y_ {t} =  beta_ {0} +  beta_ {1} t +  beta_ {2} t ^ 2 +  epsilon

Se houver duas mudanças de direção, usamos os modelos de tendência cúbica:

y_ {t} =  beta_ {0} +  beta_ {1} t +  beta_ {2} t ^ 2 +  beta_ {3} t ^ 3 +  epsilon

Por outro lado, executamos o mesmo processo que fizemos no modelo linear.

#Model variables
model_quadratic <- lm(data = df,gasoline~poly(date,2))
model_cubic <- lm(data = df,gasoline~poly(date,3))

#adjusted coefficients of determination
adj_r_squared_quadratic <- summary(model_quadratic) %>% 
  .$adj.r.squared

adj_r_squared_cubic <- summary(model_cubic) %>% 
  .$adj.r.squared

Para comparar todos os modelos, criamos uma variável que reúne todos os R ^ 2.

adj_r_squared_all <- c(linear=round(adj_r_squred_linear,4),
                   exponential=round(adj_r_squared_exponential,4),
                   quadratic=round(adj_r_squared_quadratic,4),
                   cubic=round(adj_r_squared_cubic,4))

Como podemos ver abaixo, os resultados justificam o gráfico que construímos antes; modelos polinomiais são muito melhores que modelos lineares e exponenciais. Apesar de os modelos polinomiais parecerem quase iguais, o modelo quadrático é um pouco melhor que o modelo cúbico.

adj_r_squared_all
     #linear exponential   quadratic       cubic 
     #0.6064      0.6477      0.8660      0.8653

Sazonalidade

O componente de sazonalidade representa as repetições em um período específico de tempo. Séries temporais com observações semanais mensais ou trimestrais tendem a mostrar variações sazonais que se repetem a cada ano. Por exemplo, a venda de produtos de varejo aumenta a cada ano no período de Natal ou os passeios de férias aumentam no verão. Para detectar a sazonalidade, usamos a análise de decomposição.

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br

Análise de decomposição: T, S, I, são tendência, sazonalidade e componentes aleatórios da série, respectivamente. Quando a variação sazonal aumenta à medida que a série temporal aumenta, usamos o modelo multiplicativo.

y_ {t} = T_ {t}  vezes S_ {t}  vezes I_ {t}

Se a variação parecer constante, devemos usar o modelo aditivo.

y_ {t} = T_ {t} + S_ {t} + I_ {t}

Para descobrir qual modelo é adequado, precisamos analisá-lo no gráfico.

#we create a time series variable for the plot
gasoline_ts <- ts(gasoline_df$gasoline,start = c(2013,1),end = c(2019,12),frequency = 12)
#The plot have all the components of the series(T, S, I)
plot(gasoline_ts,lwd=2,ylab="Gasoline")

Como podemos ver no gráfico acima, parece que o modelo multiplicativo seria mais adequado para os dados; especialmente se olharmos entre 2016 e 2019. Quando analisamos o gráfico, vemos algumas variações de conjuntura causada pela expansão ou contração da economia. Diferentemente da sazonalidade, as flutuações da conjuntura podem durar vários meses ou anos.

Além disso, a magnitude da flutuação da conjuntura pode mudar com o tempo, porque as flutuações diferem em magnitude e comprimento, sendo difícil refletir com dados históricos; por isso, negligenciamos o componente de conjuntura da série.

Extraindo a sazonalidade

Médias móveis (ma) costuma ser usado para separar o efeito de uma tendência da sazonalidade.

m , , termo , , movimento , , médias =  frac {Média , , de , , o , , último , , m , , observações} {m}

Usamos a média móvel cumulativa (CMA), que é a média de duas médias consecutivas, para mostrar a média móvel de ordem par. Por exemplo; dois primeiros termos ma são  bar {y} _ {6.5} e  bar {y} _ {7,5} mas não existem tais termos na série original; portanto, calculamos a média dos dois termos para encontrar o termo CMA correspondente na série:

 bar y_7 = frac { bar y_ {6,5} +  bar y_ {7,5}} {2}

O valor padrão do argumento "central" da função ma permanece TRUE para obter o valor CMA.

library(forecast)

gasoline_trend <- forecast::ma(gasoline_ts,12)

Se notado,  bar y_ {t} elimina variações sazonais e aleatórias; conseqüentemente:  bar y_ {t} = T_ {t}

Desde a y_ {t} = T_ {t}  vezes S_ {t}  vezes I_ {t} e  bar y_ {t} = T_ {t}, a variável detrend é encontrada por:  frac { bar y_ {t}} {y_ {t}} = S_ {t}  vezes I_ {t}; também é chamado método da relação entre médias móveis.

gasoline_detrend <- gasoline_ts/gasoline_trend

Cada mês possui várias proporções, onde cada proporção coincide com um ano diferente; nesta amostra, cada mês tem sete proporções diferentes. A média aritmética é usada para determinar o valor comum para cada mês; ao fazer isso, eliminamos o componente aleatório e subtraímos a sazonalidade da variável detrend. É chamado de índice sazonal não ajustado.

unadjusted_seasonality <- sapply(1:12, function(x) mean(window(gasoline_detrend,c(2013,x),c(2019,12),deltat=1),na.rm = TRUE)) %>% round(4)

Os índices sazonais para dados mensais devem ser preenchidos para 12, com uma média de 1; para fazer isso, cada índice sazonal não ajustado é multiplicado por 12 e dividido pela soma de 12 índices sazonais não ajustados.

adjusted_seasonality <- (unadjusted_seasonality*(12/sum(unadjusted_seasonality))) %>% 
                        round(4)

A média de todos os índices sazonais ajustados é 1; se um índice for igual a 1, significa que não haveria sazonalidade. Para traçar a sazonalidade:

#building a data frame to plot in ggplot
adjusted_seasonality_df <- data_frame(
  months=month.abb,
  index=adjusted_seasonality)

#converting char month names to factor sequentially
adjusted_seasonality_df$months <- factor(adjusted_seasonality_df$months,levels = month.abb)

ggplot(data = adjusted_seasonality_df,mapping = aes(x=months,y=index))+
 geom_point(aes(colour=months))+
 geom_hline(yintercept=1,linetype="dashed")+
 theme(legend.position = "none")+
 ggtitle("Seasonality of Unleaded Gasoline Prices for 95 Octane")+
 theme(plot.title = element_text(h=0.5))

Como podemos ver acima, há uma queda sazonal de aproximadamente três por cento em janeiro e dezembro; por outro lado, há um aumento sazonal de dois por cento em maio e junho. A sazonalidade não é vista em março, julho e agosto; porque seus valores de índice são aproximadamente iguais a 1.

Decomposição gráfica da série temporal

Primeiro, mostraremos a linha de tendência na série temporal.

#Trend is shown by red line
plot(gasoline_ts,lwd=2,ylab="Gasoline")+
lines(gasoline_trend,col="red",lwd=3)

E isolará a linha de tendência do gráfico:

plot(gasoline_trend,lwd=3,col="red",ylab="Trend")

Vamos ver a linha da sazonalidade:

plot(as.ts(rep(unadjusted_seasonality,12)),lwd=2,ylab="Seasonality",xlab="")

E, finalmente, mostraremos a linha da aleatoriedade; a fim de fazer isso:
I_ {t} =  frac {y_ {t}} {S_ {t}  times T_ {t}}

randomness <- gasoline_ts/(gasoline_trend*unadjusted_seasonality)
plot(randomness,ylab="Randomness",lwd=2)

Referências

Sanjiv Jaggia, Alison Kelly (2013). Inteligência de negócios: comunicação com números.

var vglnk = {key: '949efb41171ac6ec1bf7f206d57e90b8'};

(função (d, t) {
var s = d.createElement
s.src = '//cdn.viglink.com/api/vglnk.js';
var r = d.getElementsByTagName
} (documento, 'script'));

Para Deixe um comentário para o autor, siga o link e comente no blog: R-posts.com.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da dados / R.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br
Leia Também  Explorando os dados demográficos dos distritos congressionais dos EUA em R