Componente trend

Il local linear trend is that more utilizzato per the models UCM ed is 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 can essere interpretato as a retta, in the quale the coefficient angolare and l’intercetta evolvono as RW.
• Pendenza (slope): \(\beta_t = \beta_{t-1} + \zeta_t\) con \(\zeta_t \sim WN(0,\sigma^2 \zeta)\). Si tratta di a RW that evolve in the tempo and rappresenta the pendenza di a retta.

Le components are sviluppi stocastici di models deterministici. Si parte da a retta with respect a \(t = 1,2,...: \mu_t = \mu_0 +t·\beta\). Per far sì that the coefficient angolare and l’intercetta possano muoversi in the tempo it scrive the function in forma incrementale and a this it aggiunge a component rumore that influenza the retta in the tempo poiché it accumula:

\[\mu_t = \mu_{t-1}+\beta+ \eta_t \\where: \\ \beta = drift \\ \eta_t \sim WN\] We obtain così a RWD: the parte of the shock has mean 0 and can far salire or scendere the level (incremento random), mentre if \(\beta\) is positivo l’incremento sarà meanmente positivo, altrimenti sarà meanmente 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 random and pendenza fissa).
3. Se \(Var(\eta_t) = 0,Var(\zeta_t) > 0 \Rightarrow \mu_t = IRW,\) that is random walk integrato or trend liscio \(\sim I(2)\) perché the retta evolve only 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 only the 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 a boxcox but ci fidiamo and andiamo in logaritmo

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

#va very meglio, 攼㸸 a serie that pthat ismo modelszzare

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

trend<-1:n
#facciamo a before regression delthe series storica logaritmica
#considerando constant and 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:
##              Eestimatete 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 seasonal

Ci are due diversi modi per creare a component di seasonaltà in a model UCM:

• Metodo of the dummy stocastiche;
• Metodo of the sinusoidi stocastiche.

La seasonaltà is qualcosa that it ripete each dato periodo di tempo and modifica the value osservato in the singoli periodi with respect to the level medio delthe series storica, but senza influenzare tale level

#castruiamo variables factor con the 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 as sarebbe the 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:
##               Eestimatete 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

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

#costruzione sinusoide: cos(freq*tempo) & sen(freq*tempo)
ciclo<-cbind(cos(outer(1:144,frequencies_seasonal)),sin(outer(1:144,frequencies_seasonal)))
#l'ultima colonna va tolta perch攼㸸 攼㸸 always 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:
##               Eestimatete 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 can essere utilizzato as component a sè per normalizzare the ciclo economico or as parte con certe frequencies permette di normalizzare the seasonaltà.

Viene generato da a function sinusoidale ed is generato da a vector di dimensione 2 (ai fini of the previsione is indifferente prendere \(\psi_t \space or \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\) is a fattore di smorzamento.

\(R(\lambda)=\left( \begin{array}{cc} \cos \lambda & \sin \lambda \\ -\sin \lambda & \cos \lambda \end{array} \right)\) is the matrix di rotazione and \(\lambda\) is the frquenza of the ciclo stocastico (angolo of the sinusoide).

\(\left( \begin{array}{cc} k_t \\ k_t^* \end{array} \right)\) are the rispettivi WN di \(\psi_t\) and di \(\psi_t^*\)

Per arrivare a this formula finale partiamo from the definition di a sinusoide that varia in the tempo con ampiezza R frequency \(\lambda\) and con fase \(\phi\):

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

Quindi the sinusoide dipende da 3 fattori:

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

La frequency indica quanto velocemente the function ci mette a fare a giro completo, quindi more \(\lambda\) is grande maggiore sarà the numero di oscillazioni in a unità di tempo. Per definire in modo more chiaro the frequency definiamo \(\lambda\) as \(\lambda = \frac{2 \space \pi}{periodo}\) and the periodo is the numero di istanti temporal that ci mette the sinusoide per completare the suo ciclo.

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

R modifica the range di oscillazione (ampiezza) of the sinusoide:

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

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

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

This forma funzionale is very utile per capire l’utilità di each singolo parameter but given that \(\phi \space and \space \lambda\) are all’interno of the coseno in the calculation of the regressori sarebbe necessary minimizzare the quadrati in modo not linear, per this spesso it utilizza un’altra scrittura of the sinusoide that is 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 this modo considerando \(\cos(\lambda \space t) \space and \space \sin ( \lambda \space t)\) as due regressori A and B it possono estimatere con the minimi quadrati ordinari. Faciendo the Stesso ragionamento fatto sopra con the coseno per the seno of the sinusoide troviamo the matrix di rotazione.

La matrix, quindi, \(R(\lambda)\) is a matrix quadrata ortogonale di rango 2 that permette the rotazione in senso orario in a piano di 2 dimensioni di a angolo \(\lambda\) intorno all’origine.

Infine andiamo ad analizzare l’utilità of the parameter presente in the formula iniziale \(\rho\). Andiamo ad analizzare the situazioni possibili and the possibili risultati:

• \(\rho=0 \Rightarrow\) implica that \(\psi_t \sim WN\);
• \(\rho=1 \Rightarrow\) ciclo not stationary;
• \(0 < \rho<1 \Rightarrow\) the ciclo economico has a oscillazione di medio-breve termine (solitamente situazione more plausibile) and \(\rho\) serve per rendere the ciclo stationary dove the 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 a example con a model strutturato as 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 the chiameremo rho pewrche R chiamiamo the function 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 the function di rotazione
R = function ( lambda ){
  co = cos(lambda)
  it = 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:
# example 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 as sarebbe the 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:
##               Eestimatete 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 the time series presentano cambiamenti repentini and if ne possono considerare quattro tipi:

1. Outlier additivo (AO): in a certa data the series storica assume a value very diverso dagli altri but from the data successiva torna a seguire l’andamento that aveva precedentemente (ad es. a evento that dura only a giorno if the series is 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 a serie di date successive (potrebbe essere dovuto ad a evento that dura more di a giorno but poi it 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 a cambio di level definitivo (variable 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 delthe series 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}\]

This model presuppone that \(X\) to the tempo \(t\) abbia effetto only su \(Y\) to the tempo \(t\) and that the suo effetto both nullo in the tempi successivi. Si tratta di a situazione riduttiva perché ad es. a inveestimatesnto in pubblicità has both effetti immeanti that duraturi, also if it affievoliscono in the tempo.

Regressori dinamici

Il model con the regressori statici presuppone that \(X\) to the tempo \(t\) abbia effetto only su \(Y\) to the tempo \(t\) and that the suo effetto both nullo in the tempi successivi. Si tratta di a situazione riduttiva perché ad es. a inveestimatesnto in pubblicità has both effetti immeanti that duraturi, also if it affievoliscono in the tempo. Potrebbe quindi essere utile introdurre of the lags.