Componente trend

Il local linear trend è quello più utilizzato per i modelli UCM ed è formmato da due componemti:

• Livello: \(\mu_t = \mu_{t-1} + \beta_{t-1} + \eta_t\) con \(\eta_t \sim WN(0,\sigma^2 \eta)\). Localmente può essere interpretato come una retta, nella quale il coefficiente angolare e l’intercetta evolvono come RW.
• Pendenza (slope): \(\beta_t = \beta_{t-1} + \zeta_t\) con \(\zeta_t \sim WN(0,\sigma^2 \zeta)\). Si tratta di un RW che evolve nel tempo e rappresenta la pendenza di una retta.

Le componenti sono sviluppi stocastici di modelli deterministici. Si parte da una retta rispetto a \(t = 1,2,...: \mu_t = \mu_0 +t·\beta\). Per far sì che il coefficiente angolare e l’intercetta possano muoversi nel tempo si scrive la funzione in forma incrementale e a questa si aggiunge una componente rumore che influenza la retta nel tempo poiché si accumula:

\[\mu_t = \mu_{t-1}+\beta+ \eta_t \\dove: \\ \beta = drift \\ \eta_t \sim WN\] Si ottiene così un RWD: la parte degli shock ha media 0 e può far salire o scendere il livello (incremento casuale), mentre se \(\beta\) è positivo l’incremento sarà mediamente positivo, altrimenti sarà mediamente negativo (incremento deterministico).

1. Se \(Var(\eta_t) = Var(\zeta_t) = 0\Rightarrow \mu_t = retta\sim I(0)\).
2. Se \(Var(\eta_t) > 0,Var(\zeta_t) = 0\Rightarrow \mu_t = RWD\sim I(1)\) (intercetta casuale e pendenza fissa).
3. Se \(Var(\eta_t) = 0,Var(\zeta_t) > 0 \Rightarrow \mu_t = IRW,\) cioè random walk integrato o trend liscio \(\sim I(2)\) perché la retta evolve solo tramite cambiamenti di pendenza.
4. Se \(Var(\eta_t) > 0,Var(\zeta_t) = 0,\beta0 = 0\Rightarrow \mu_t = RW\sim I(1)\) (local level model).

#estraiamo l'anno
anno<-floor(time(AirPassengers)) #estraiamo solo la parte intera

medie<-tapply(AirPassengers, anno, mean); medie

##     1949     1950     1951     1952     1953     1954     1955     1956 
## 126.6667 139.6667 170.1667 197.0000 225.0000 238.9167 284.0000 328.2500 
##     1957     1958     1959     1960 
## 368.4167 381.0000 428.3333 476.1667

devstandard<-tapply(AirPassengers, anno, sd); devstandard

##     1949     1950     1951     1952     1953     1954     1955     1956 
## 13.72015 19.07084 18.43827 22.96638 28.46689 34.92449 42.14046 47.86178 
##     1957     1958     1959     1960 
## 57.89090 64.53047 69.83010 77.73713

plot(medie, devstandard)
abline(lm(devstandard~medie), col=2)

#bisognerebbe applicare un boxcox ma ci fidiamo e andiamo in logaritmo

y<-log(AirPassengers)
plot.ts(y)

#va molto meglio, 攼㸸 una serie che possiamo modellizzare

#partiamo con il trend:
n<-length(y)

trend<-1:n
#facciamo una prima regressione della serie storica logaritmica
#considerando costante e trend

reg1<-lm(y~trend)
plot(as.numeric(y), type="l")
abline(reg1, col="red")

summary(reg1)

## 
## Call:
## lm(formula = y ~ trend)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30858 -0.10388 -0.01796  0.09738  0.29538 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.8136683  0.0232940  206.65   <2e-16 ***
## trend       0.0100484  0.0002787   36.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.139 on 142 degrees of freedom
## Multiple R-squared:  0.9015, Adjusted R-squared:  0.9008 
## F-statistic:  1300 on 1 and 142 DF,  p-value: < 2.2e-16

Componente stagionale

Ci sono due diversi modi per creare una componente di stagionalità in un modello UCM:

• Metodo delle dummy stocastiche;
• Metodo delle sinusoidi stocastiche.

La stagionalità è qualcosa che si ripete ogni dato periodo di tempo e modifica il valore osservato nei singoli periodi rispetto al livello medio della serie storica, ma senza influenzare tale livello

#castruiamo variabili factor con i nostri mesi:
mesi<-c("gen","feb","mar","apr","magg","giugno","lugl","ago","sett","ott","nov","dec")
mesi<-factor(rep(mesi,(1960-1949)+1),levels = c("gen","feb","mar","apr","magg","giugno","lugl","ago","sett","ott","nov","dec"))

head(model.matrix(y~mesi))

##   (Intercept) mesifeb mesimar mesiapr mesimagg mesigiugno mesilugl mesiago
## 1           1       0       0       0        0          0        0       0
## 2           1       1       0       0        0          0        0       0
## 3           1       0       1       0        0          0        0       0
## 4           1       0       0       1        0          0        0       0
## 5           1       0       0       0        1          0        0       0
## 6           1       0       0       0        0          1        0       0
##   mesisett mesiott mesinov mesidec
## 1        0       0       0       0
## 2        0       0       0       0
## 3        0       0       0       0
## 4        0       0       0       0
## 5        0       0       0       0
## 6        0       0       0       0

reg2<-lm(y~trend+mesi)


#vediamo come sarebbe il nostro fit:
plot(as.numeric(y), type="l")
lines(reg2$fitted.values, col="red")

summary(reg2)

## 
## Call:
## lm(formula = y ~ trend + mesi)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.156370 -0.041016  0.003677  0.044069  0.132324 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7267804  0.0188935 250.180  < 2e-16 ***
## trend        0.0100688  0.0001193  84.399  < 2e-16 ***
## mesifeb     -0.0220548  0.0242109  -0.911  0.36400    
## mesimar      0.1081723  0.0242118   4.468 1.69e-05 ***
## mesiapr      0.0769034  0.0242132   3.176  0.00186 ** 
## mesimagg     0.0745308  0.0242153   3.078  0.00254 ** 
## mesigiugno   0.1966770  0.0242179   8.121 2.98e-13 ***
## mesilugl     0.3006193  0.0242212  12.411  < 2e-16 ***
## mesiago      0.2913245  0.0242250  12.026  < 2e-16 ***
## mesisett     0.1466899  0.0242294   6.054 1.39e-08 ***
## mesiott      0.0085316  0.0242344   0.352  0.72537    
## mesinov     -0.1351861  0.0242400  -5.577 1.34e-07 ***
## mesidec     -0.0213211  0.0242461  -0.879  0.38082    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0593 on 131 degrees of freedom
## Multiple R-squared:  0.9835, Adjusted R-squared:  0.982 
## F-statistic: 649.4 on 12 and 131 DF,  p-value: < 2.2e-16

frequenze_stagionali<-2*pi*(1:6)/12

#costruzione sinusoide: cos(freq*tempo) & sen(freq*tempo)
ciclo<-cbind(cos(outer(1:144,frequenze_stagionali)),sin(outer(1:144,frequenze_stagionali)))
#l'ultima colonna va tolta perch攼㸸 攼㸸 sempre 0

reg3<-lm(y~trend+ciclo)


plot(as.numeric(y), type="l")
lines(reg3$fitted.values, col="red")

summary(reg3)

## 
## Call:
## lm(formula = y ~ trend + ciclo)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.157681 -0.041442  0.003145  0.044222  0.131433 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.812e+00  1.000e-02 481.015  < 2e-16 ***
## trend        1.007e-02  1.198e-04  84.043  < 2e-16 ***
## ciclo1      -1.418e-01  7.027e-03 -20.174  < 2e-16 ***
## ciclo2      -2.272e-02  7.049e-03  -3.223 0.001603 ** 
## ciclo3       2.730e-02  7.017e-03   3.891 0.000158 ***
## ciclo4       2.212e-02  7.020e-03   3.150 0.002024 ** 
## ciclo5       5.535e-03  7.018e-03   0.789 0.431702    
## ciclo6       3.423e-03  6.094e-03   0.562 0.575296    
## ciclo7      -4.936e-02  7.030e-03  -7.022 1.09e-10 ***
## ciclo8       7.863e-02  7.026e-03  11.192  < 2e-16 ***
## ciclo9      -8.823e-03  7.048e-03  -1.252 0.212892    
## ciclo10      2.560e-02  7.016e-03   3.649 0.000380 ***
## ciclo11      2.138e-02  7.016e-03   3.047 0.002797 ** 
## ciclo12      5.377e+10  3.909e+11   0.138 0.890801    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05953 on 130 degrees of freedom
## Multiple R-squared:  0.9835, Adjusted R-squared:  0.9818 
## F-statistic:   595 on 13 and 130 DF,  p-value: < 2.2e-16

Componente ciclica

Il ciclo stocastico può essere utilizzato come componente a sè per normalizzare il ciclo economico o come parte con certe frequenze permette di normalizzare la stagionalità.

Viene generato da una funzione sinusoidale ed è generato da un vettore di dimensione 2 (ai fini della previsione è indifferente prendere \(\psi_t \space o \space \psi_t^*\):

\[\left( \begin{array}{cc} \psi_t \\ \psi_t^* \end{array} \right) = \rho \space \left( \begin{array}{cc} \cos \lambda & \sin \lambda \\ -\sin \lambda & \cos \lambda \end{array} \right) \left( \begin{array}{cc}\psi_{t-1} \\ \psi_{t-1}^* \end{array} \right) + \left( \begin{array}{cc} k_t \\ k_t^* \end{array} \right) \] Dove:

\(\rho\) è un fattore di smorzamento.

\(R(\lambda)=\left( \begin{array}{cc} \cos \lambda & \sin \lambda \\ -\sin \lambda & \cos \lambda \end{array} \right)\) è la matrice di rotazione e \(\lambda\) è la frquenza del ciclo stocastico (angolo della sinusoide).

\(\left( \begin{array}{cc} k_t \\ k_t^* \end{array} \right)\) sono i rispettivi WN di \(\psi_t\) e di \(\psi_t^*\)

Per arrivare a questa formula finale partiamo dalla definizione di una sinusoide che varia nel tempo con ampiezza R frequenza \(\lambda\) e con fase \(\phi\):

\[ f(t) = R \cos (\phi + \lambda \space t) \]

Quindi la sinusoide dipende da 3 fattori:

• ampiezza \(R\)
• frequenza $$
• fase \(\phi\)

La frequenza indica quanto velocemente la funzione ci mette a fare un giro completo, quindi più \(\lambda\) è grande maggiore sarà il numero di oscillazioni in una unità di tempo. Per definire in modo più chiaro la frequenza definiamo \(\lambda\) come \(\lambda = \frac{2 \space \pi}{periodo}\) e il periodo è il numero di istanti temporali che ci mette la sinusoide per completare il suo ciclo.

Vediamo alcuni esempi (per adesso ipotiziamo \(R=1 \space e \space \phi=0\):

R modifica il range di oscillazione (ampiezza) della sinusoide:

Vediamo alcuni esempi (per adesso ipotiziamo \(\lambda=\pi \space e \space \phi=0\):

Infine la fase modifica il punto da cui parte il coseno, quindi sposta il grafico di una quantità pari a \(-\phi/\lambda\).

Vediamo alcuni esempi (per adesso ipotiziamo \(R=1 \space e \space \lambda=\pi\):

Questa forma funzionale è molto utile per capire l’utilità di ogni singolo parametro ma dato che \(\phi \space e \space \lambda\) sono all’interno del coseno nel calcolo dei regressori sarebbe necessario minimizzare i quadrati in modo non lineare, per questo spesso si utilizza un’altra scrittura della sinusoide che è equivalente:

\[ f(t) = R \cos(\phi + \lambda \space t) = \\ = R \cos(\phi) \cos(\lambda \space t) - R \sin(\phi) \sin(\lambda \space t) = \\ = A \cos(\lambda \space t) + B \sin ( \lambda \space t) \]

In questo modo considerando \(\cos(\lambda \space t) \space e \space \sin ( \lambda \space t)\) come due regressori A e B si possono stimare con i minimi quadrati ordinari. Faciendo lo Stesso ragionamento fatto sopra con il coseno per il seno della sinusoide troviamo la matrice di rotazione.

La matrice, quindi, \(R(\lambda)\) è una matrice quadrata ortogonale di rango 2 che permette la rotazione in senso orario in un piano di 2 dimensioni di un angolo \(\lambda\) intorno all’origine.

Infine andiamo ad analizzare l’utilità del parametro presente nella formula iniziale \(\rho\). Andiamo ad analizzare le situazioni possibili e i possibili risultati:

• \(\rho=0 \Rightarrow\) implica che \(\psi_t \sim WN\);
• \(\rho=1 \Rightarrow\) ciclo non stazionario;
• \(0 < \rho<1 \Rightarrow\) il ciclo economico ha una oscillazione di medio-breve termine (solitamente situazione più plausibile) e \(\rho\) serve per rendere il ciclo stazionario dove i primi due modenti sono:\[ E(\psi_t)=0\] \[ Var(\psi_t)=E(\psi_t \psi_t^T) = \frac{\sigma_k^2}{1- \rho^2} I_2\]

Andiamo a fare un esempio con un modello strutturato come segue:

\[Y_t= \mu_t+\psi_t + \epsilon_t \\.\\ Dove: \\ \mu_t = \mu_{t-1}+\beta+ \eta_t \space \space con \space \space \eta_t \sim WN(0, \sigma^2_{\eta})\\ \left( \begin{array}{cc} \psi_t \\ \psi_t^* \end{array} \right) = \rho \space \left( \begin{array}{cc} \cos \lambda & \sin \lambda \\ -\sin \lambda & \cos \lambda \end{array} \right) \left( \begin{array}{cc}\psi_{t-1} \\ \psi_{t-1}^* \end{array} \right) + \left( \begin{array}{cc} k_t \\ k_t^* \end{array} \right) \space \space con \space \space \left( \begin{array}{cc} k_t \\ k_t^* \end{array} \right) \sim WN(\underline{0}, \sigma_k^2 \space I_2) \\ \epsilon_t = WN(0, \sigma^2_{\epsilon})\]

y<-log(AirPassengers)
trend<-1:144

R=1 # da adesso lo chiameremo rho pewrche R chiamiamo la funzione di rotazione
rho=1
phi=0
periodo=12
lambda=2*pi/periodo
x=seq(1,length(AirPassengers),1)
ciclo<- matrix(c(0,1),2,144)

#Creiamo la funzione di rotazione
R = function ( lambda ){
  co = cos(lambda)
  si = sin(lambda)
  matrix(c(co,-si,si,co),2,2)
}


sig_kappa=0.01
for ( t in 2:n){
  ciclo[,t]=rho*R(lambda) %*% ciclo[,t-1]+rnorm(2,sd=sig_kappa)
}



#per cicli di ordine superiore usare:
# esempio ciclo di ordine 2
#for ( t in 2:n){
#  ciclo1[,t]=rho*R(lambda) %*% ciclo1[,t-1]+rnorm(2,sd=sig_kappa)
#  ciclo2[,t]=rho*R(lambda) %*% ciclo2[,t-1]+ciclo1[,t-1]  #ciclo di ordine superiore
#}


reg4<-lm(y~trend+ciclo[1,]+ciclo[2,])


#vediamo come sarebbe il nostro fit:
plot(as.numeric(y), type="l")
lines(reg4$fitted.values, col="red")

summary(reg4)

## 
## Call:
## lm(formula = y ~ trend + ciclo[1, ] + ciclo[2, ])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.261856 -0.065582  0.001491  0.068906  0.177989 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.8147793  0.0150975 318.911  < 2e-16 ***
## trend        0.0100303  0.0001807  55.493  < 2e-16 ***
## ciclo[1, ]   0.0350801  0.0102077   3.437 0.000776 ***
## ciclo[2, ]  -0.1392321  0.0101752 -13.683  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08996 on 140 degrees of freedom
## Multiple R-squared:  0.9593, Adjusted R-squared:  0.9585 
## F-statistic:  1101 on 3 and 140 DF,  p-value: < 2.2e-16

Regressori statici

Spesso le serie storiche presentano cambiamenti repentini e se ne possono considerare quattro tipi:

1. Outlier additivo (AO): in una certa data la serie storica assume un valore molto diverso dagli altri ma dalla data successiva torna a seguire l’andamento che aveva precedentemente (ad es. un evento che dura solo un giorno se la serie è giornaliera). \[AO_t = \begin{cases} 1, & \mbox{if } \space t=t_0 \\ 0, & \mbox{if } \space t \ne t_0 \end{cases}\]

2. Temporary change (TC): avviene su una serie di date successive (potrebbe essere dovuto ad un evento che dura più di un giorno ma poi si esaurisce).\[TC_t = \begin{cases} 1, & \mbox{if } \space t_0 \le t \le t_1 \\ 0, & \space altrimenti \end{cases} \]

3. Level shift (LS): avviene un cambio di livello definitivo (variabile scalino)\[LS_t = \begin{cases} 1, & \mbox{if } \space t\ge t_0 \\ 0, & \mbox{if } \space t < t_0 \end{cases}\]
4. Slope shift (SS): cambio di pendenza permanente della serie stori\[SS_t = \begin{cases} t-t_0+1, & \mbox{if } \space t\ge t_0 \\ 0, & \mbox{if } \space t < t_0 \end{cases}\]

Questo modello presuppone che \(X\) al tempo \(t\) abbia effetto solo su \(Y\) al tempo \(t\) e che il suo effetto sia nullo nei tempi successivi. Si tratta di una situazione riduttiva perché ad es. un investimento in pubblicità ha sia effetti immediati che duraturi, anche se si affievoliscono nel tempo.

Regressori dinamici

Il modello con i regressori statici presuppone che \(X\) al tempo \(t\) abbia effetto solo su \(Y\) al tempo \(t\) e che il suo effetto sia nullo nei tempi successivi. Si tratta di una situazione riduttiva perché ad es. un investimento in pubblicità ha sia effetti immediati che duraturi, anche se si affievoliscono nel tempo. Potrebbe quindi essere utile introdurre dei ritardi.