Modello lineare multivariato

Vogliamo stimare il seguente modello lineare:

\[Y_t = \beta_0 + \beta_1 X_{t. 1} + \beta_2 X_{t. 2} + \epsilon_t\] Notiamo che:

• è stata introdotta l’intercetta \(\beta_0\) che rappresenta la media non condizionata di \(Y\), per costruzione (avendo simulato \(Y\) come una normale con media zero) avremo che \(\mathbb{E}\{Y\} = \beta_0 = 0\). Una volta stimato il modello ci aspettiamo che \(\beta_0 \approx 0\) e che il relativo test \(t\) non sia significativo (e.g. non possiamo rifiutare l’ipotesi nulla \(H_{0}: \beta_0 = 0\)).

• Come è rappresentata invece la media condizionata di Y dati \(X_1\) e \(X_2\)? Assumendo che la forma funzionale di Y sia quella espressa dal modello, possiamo scriverla come:

\[\mathbb{E}\{Y | X_1, X_2\} = \beta_0 + \beta_1 X_{1} + \beta_2 X_{2}\].

• Dove è finito \(\epsilon\)? Per costruzione \(\mathbb{E}\{\epsilon | X_1, X_2\} = 0\)

model <- lm(Y ~ X1 + X2, data = as.data.frame(X))
summary(model)

## 
## Call:
## lm(formula = Y ~ X1 + X2, data = as.data.frame(X))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.68049 -0.65667 -0.02984  0.68851  2.57767 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.007202   0.030605   0.235    0.814    
## X1           0.302251   0.029943  10.094  < 2e-16 ***
## X2          -0.214021   0.029245  -7.318 5.17e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9674 on 997 degrees of freedom
## Multiple R-squared:  0.1319, Adjusted R-squared:  0.1301 
## F-statistic: 75.74 on 2 and 997 DF,  p-value: < 2.2e-16

x <- seq(-3, 3, 0.01)
p <- stats::dt(x, df = 999)
ggplot()+
  geom_line(aes(x, p)) +
  geom_segment(aes(x = 0.235, xend=0.235, y=0, yend=stats::dt(0.235, df = 999)), col="red")+
  geom_segment(aes(x = -0.235, xend=-0.235, y=0, yend=stats::dt(-0.235, df = 999)), col="red")+
  geom_label(aes(x = 1, y = 0.08, label = "P(t value > t)"))+
  geom_label(aes(x = -1, y = 0.08, label = "P(t value < -t)"))+
  scale_x_continuous(breaks = c(seq(-3, -1, 1),-0.235, 0.235, seq(1, 3, 1)))+
  xlab("t value")+
  ylab("probability")+
  ggtitle("Test t sull'intercetta")

# coda di destra: P(t value > t)
p_dx <- stats::pt(0.235321, df = 999) - 0.5
# coda di sinistra: P(t value < -t)
p_sx <- 0.5 - stats::pt(-0.235321, df = 999) 
# area centrale: 1 - coda sx - coda dx = P(t value > |t|)
1 - (p_sx + p_dx)

## [1] 0.8140078

model <- lm(Y ~ X1 + X2 - 1, data = as.data.frame(X))
summary(model)

## 
## Call:
## lm(formula = Y ~ X1 + X2 - 1, data = as.data.frame(X))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.67297 -0.64963 -0.02289  0.69565  2.58492 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## X1  0.30231    0.02993  10.101  < 2e-16 ***
## X2 -0.21423    0.02922  -7.332 4.68e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9669 on 998 degrees of freedom
## Multiple R-squared:  0.132,  Adjusted R-squared:  0.1303 
## F-statistic: 75.92 on 2 and 998 DF,  p-value: < 2.2e-16

# valori fittati
Y_hat <- predict(model)

# residui
epsilon <- residuals(model)

Test di omoschedasticità: Test di White

Per effettuare il test di white fittiamo un’ulteriore regressione dei quadrati dei residui sui quadrati delle variabili indipendenti (\(X_1\) e \(X_2\)):

\[\hat{\epsilon}^2_t = \gamma_0 + \gamma_{1} X_{1,t}^2 + \gamma_{2} X_{2,t}^2 + \gamma_{3} X_{1,t} X_{2,t}\]

• L’intercetta \(\gamma_0\) deve essere inclusa

# quadrati dei residui
epsilon_2 = epsilon^2
# quadrato di X1
x1_2 = X[,2]^2
# quadrato di X2
x2_2 = X[,3]^2
# prodotto incrociato di X1 e X2
x1_x2 = X[,2] * X[,3]
# regressione per il test di white
lm(epsilon_2 ~ x1_2 + x2_2 + x1_x2) %>%
  summary()

## 
## Call:
## lm(formula = epsilon_2 ~ x1_2 + x2_2 + x1_x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9947 -0.8347 -0.4722  0.3665  6.2107 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.948004   0.055919  16.953   <2e-16 ***
## x1_2        -0.018032   0.026347  -0.684    0.494    
## x2_2         0.004039   0.025247   0.160    0.873    
## x1_x2       -0.019325   0.038330  -0.504    0.614    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.247 on 996 degrees of freedom
## Multiple R-squared:  0.0008076,  Adjusted R-squared:  -0.002202 
## F-statistic: 0.2683 on 3 and 996 DF,  p-value: 0.8482

Ora dobbiamo osservare la statistica F (F-statistic), che è 0.2683 e il relativo p-value 0.8482. Possiamo affermare che non rifiutiamo l’ipotesi nulla, che rappresenta la situazione di omoschedasticità, con ogni livello di significatività (> 10%). Ricordiamo che la statistica F riportata saggia l’ipotesi nulla:

\[H_0: \gamma_1 = \gamma_2 = \gamma_3 = 0 \; \; \; \text{or} \; \; \; H_0: \sigma_1^2 = \sigma_2^2 = \sigma_3^2 = \sigma^2 = \gamma_0 \approx \text{homoskedasticity}\] ovvero siamo in presenza di omoschedasticità.

Test di omoschedasticità: Breush-Pagan

Vediamo ora se anche il test di Breush-Pagan ci fa arrivare alla stessa conclusione:

lmtest::bptest(model)

## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 3.0367, df = 1, p-value = 0.0814

Come possiamo vedere il p-value assoociato è 0.0814, di conseguenza non rifiutiamo l’ipotesi nulla \(H_0\) (omoschedasticità) con il 5% di significatività.

Test di incorrelazione seriale: Durbin-Watson

L’idea del test di Durbin-Watson W è di stimare il seguente modello sui residui \(\hat{\epsilon_t}\):

\[\hat{\epsilon}_t = \alpha \hat{\epsilon}_{t-1} + u_t \; \; \; , \; \; \; u_t \sim i.i.d.(0, \sigma^2) \] La statistica test, spesso indicata con DW, si calcola:

\[\text{DW} = \frac{\sum_{t=2}^{n} ( \hat{\epsilon}_{t} - \hat{\epsilon}_{t-1})^{2} }{\sum_{t=2}^{n} \hat{\epsilon}_{t}^2} \approx 2(1 - \alpha) \]

epsilon_lag_1 = lag(epsilon)
# statistica test di durbin-watson 
dw <- sum((epsilon-epsilon_lag_1)^2, na.rm = TRUE)/(sum(epsilon[-1]^2))
data.frame(dw = dw) %>%
  knitr::kable() %>%
  kableExtra::kable_classic_2() %>%
  kableExtra::kable_styling(latex_options = "hold_position")

Sotto l’ipotesi nulla \(H_0: \alpha = 0\) avremo che \(DW \approx 2(1-0) = 2\). Il test genera una statistica test sempre compresa tra 0 e 4, tuttavia pur calcolando la statistica test, per stabilire se rifiutare o meno l’ipotesi nulla occorre conoscere i valori limite, esprimibili attraverso tavole statistiche. Essendo il loro calcolo non molto intuitivo spesso i software riportano il p-value gia calcolato.

lmtest::dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 2.0888, p-value = 0.9247
## alternative hypothesis: true autocorrelation is greater than 0

Essendo il p-value 0.9247 non rifiutiamo l’ipotesi nulla, quindi siamo in presenza di assenza di autocorrelazione.

Test di incorrelazione seriale: Test Pormanteau

Il Test Pormanteau si basa sull’assunzione che i residui del modello siano \(\epsilon_t \sim i.i.d.(0, \sigma^2)\), ovvero la varianza dei residui è costante (omoschedasticità). Definiamo il coefficiente di autocorrelazione campionario al ritardo \(k\) come:

\[\hat{\rho}_k = \frac{\sum_{t=2}^{k} \hat{\epsilon}_{t} \hat{\epsilon}_{t-k}}{\sum_{t=k}^{n} \hat{\epsilon}_{t}^2}\] Possiamo costruire il vettore contenente \(m\)-correlazioni che converge in distribuzione a una normale multivariata:

\[\sqrt{n} \begin{bmatrix} \hat{\rho}_1 \\ \hat{\rho}_2 \\ \vdots \\ \hat{\rho}_m \end{bmatrix} \overset{d}{\rightarrow} \mathcal{N}(\underset{m \times 1}{0} \; , \; \underset{m\times m}{ \mathbb{I}})\] Da questo risultato, sapendo che la somma dei quadrati di \(m\) variabili normali è distribuita come un \(\chi(m)\) possiamo costruire la statistica test \(Q_m\):

\[Q_m = \sum_{k = 1}^{m} \hat{\rho}_{k}^{2} \overset{d}{\rightarrow} \chi(m)\]

test_pormanteau <- function(x, m = 1){
  
  n <- length(x) # numero osservazioni 

  # funzione per il calcolo della correlazione
  rho <- function(x, k = 1){

    # lag k
    x_lag_k = lag(x, k)
    
    # correlation 
    sum(x*x_lag_k, na.rm = TRUE)/sum(x[-c(1:k)]^2)
    }
  
  # calcolo delle correlazioni da 1 a k 
  vec_rho <- c()
  for(i in 1:m){
    vec_rho[i] <-  rho(x, k = i)
  }
  # statistica test 
  Qm <- n*sum(c(vec_rho)^2)
  # pvalue 
  pvalue <- 1 - pchisq(Qm, df = m)
  
  data.frame(Qm = Qm, m = m, pvalue = pvalue)
}

test_pormanteau(epsilon, 5) %>%
  knitr::kable() %>%
  kableExtra::kable_classic_2() %>%
  kableExtra::kable_styling(latex_options = "hold_position")

In questo caso abbiamo effettuato il test per \(m = 5\) ritardi ottenendo un p-value di 0.685588. Di conseguenza non rifiutiamo l’ipotesi nulla con il 5% di significatività (p-value = 0.068 > 0.05), tuttavia la rifiutiamo all’10% di significatività (p-value = 0.068 < 0.10).