## Simuliamo violazioni delle ipotesi OLS

## La variabile x sarà comune a tutti gli scenari simulati, 
## e così a e b
x <- rnorm(100)
a <- 10
b <- 0.6

## scriviamo una funzione "comoda" per visualizzare
plotres <- function(mod, ...) {
  plot(resid(mod), col="red", type="l", ...)
  abline(h=0)
}

## Normalità (v. Fig. 6.8)

## Caso "giusto"
e <- rnorm(100)
y <- a + b*x + e
plotres(lm(y~x))

## Caso di violazione dell'ipotesi

## a) code "grosse"
e <- rt(100, df=4)
y <- a + b*x + e
plotres(lm(y~x))
qqnorm(resid(lm(y~x)))

## b) asimmetria
e <- rchisq(100, df=4)
y <- a + b*x + e
plotres(lm(y~x))
qqnorm(resid(lm(y~x)))



## Omoschedasticità (v. Fig. 6.9)

## residui eteroschedastici in due gruppi:
e1 <- rnorm(70, sd=0.5)
e2 <- rnorm(30, sd=2)
e <- c(e1, e2)

y <- a + b*x + e
plotres(lm(y~x))


## Autocorrelazione (v. Fig. 6.11)

## positiva:
e <- rep(NA, 100)
e[1] <- rnorm(1)
for(i in 2:100) e[i] <- 0.7*e[i-1] + rnorm(1)

y <- a + b*x + e
plotres(lm(y~x))

## negativa:
e <- rep(NA, 100)
e[1] <- rnorm(1)
for(i in 2:100) e[i] <- -0.7*e[i-1] + rnorm(1)

y <- a + b*x + e
plotres(lm(y~x))


## Simuliamo errate specificazioni del modello lineare:

## Modello senza intercetta (v. Fig. 6.12)
e <- rnorm(100)
y <- a + b*x + e
plotres(lm(y~x-1), ylim=c(0,13))


## Esplicativa omessa (v. Fig. 6.13)

## a) incorrelata con x
e <- rnorm(100)
z <- rnorm(100)
y <- a + b*x + 3*z + e
plotres(lm(y~x))

## b) correlata con x
z <- runif(100, min=0, max=1)*x + rnorm(100) 
y <- a + b*x + 3*z + e
plotres(lm(y~x))

## c) trending
z <- runif(100, min=0, max=1)*1:100/10 + rnorm(100) 
y <- a + b*x + 3*z + e
plotres(lm(y~x))


## Cambiamenti strutturali (v. Fig. 6.14)
a1 <- 10
a2 <- c(rep(-5, 70), rep(0,30))
y <- (a1 + a2) + b*x + e
plotres(lm(y~x))


## Relazioni non lineari

## a) quadratica (a<0)
y <- a + b*x -0.5*x^2 + e
plotres(lm(y~x))

## b) esponenziale
y <- a + b*exp(x) + e
plotres(lm(y~x))
plot(residuals(lm(y~x))[order(x)])