Solitamente alcuni dei parametri dei nostri modelli in forma state space sono ignoti e devono essere stimati usando i dati. Se l’ipotesi di gaussianità dei disturbi del sistema è plausibile allora il filtro di Kalman consente di calcolare la verosimiglianza per qualunque set di valori dei parametri da stimare e, pertanto, è possibile stimare i parametri ignoti massimizzando numericamente la (log) verosimiglianza rispetto ai parametri. Se l’ipotesi di gaussianità invece non è sostenuta dai dati, allora le stime di massima verosimiglianza sono impossibili o comunque molto complesse da ottenere. Tuttavia massimizzare la (log) verosimiglianza gaussiana anche quando i dati non sono normali porta a delle stime il cui comportamento asintotico è comunque noto (sotto condizioni di regolarità). Tali stime prendono il nome di stime di quasi massima verosimiglianza (QML) o stime di pseudo massima verosimiglianza.
In generale, quando si devono produrre delle stime di massima verosimiglianza numerica si deve procedere nel seguente modo:
1. si scrive una funzione obiettivo, per esempio objective(pars, data) che, per ogni valore del vettore dei parametri da stimare pars e noti i dati data, restituisce il valore della log-verosimiglianza;
2. si assegnano dei valori iniziali ai parametri ignoti e si passano insieme alla funzione obiettivo a una funzione ottimizzatrice (in R solitamente si usa la funzione optim(), ma ce ne sono anche altre) che usando algoritmi di vario genere (per i nostri problemi si usano algoritmi di tipo Newton o quasi-Newton, spesso si adopera un algoritmo di nome BFGS) per aggiornare i parametri da stimare fino a quando un massimo (approssimativo) della funzione obiettivo viene raggiunto.
E’ importante ricordare che sotto condizioni di regolarità, il comportamento asintotico delle stime di massima verosimiglianza è noto (per comportamento asintotico usare: Metodo Delta e Informazione di Fisher).
Per la stime ci sono due alternative. La prima, più artigianale nell’implementazione, consente stime computazionalmente più efficienti. La seconda invece è più comoda, ma leggermente più lenta nell’esecuzione.
(Onestamente non so perchè differiscono [la differenza è minima] dovrebbero essere esattamente identiche io preferisco la seconda per la comodità e per la velocità di implementazione detto questo in alcune situazioni come quando presente un ciclo stocastico è comunque necessario utilizzare quella più artigianale)
Per i nostri fini l’unico argomento formale della funzione logLik(object) è un oggetto di classe SSModel, cioè la rappresentazione KFAS di un modello in forma state space. Tale funzione restituisce la log-verosimiglianza del modello in forma state space. Ovviamente, se alcuni parametri del modello sono da stimare dobbiamo creare una funzione obiettivo che prende i parametri da stimare e restituisce il valore della log-verosimiglianza. Tale funzione deve poi essere massimizzata numericamente.
Prendiamo come esempio la serie storica Nile e andiamo a stimare un modello con random walk più rumore (N.B. leggere con attenzione il codice ed i commenti contenuti all’interno).
Iniziamo caricando la libreria e i dati e creando il modello RW + WN per il Nile:
library(KFAS)
data(Nile)
# Modello in forma state space
modello1 <- SSModel(Nile~0+SSMtrend(1, NA), H = NA)Creo una funzione che stima le tre varianza del modello LLT + WN :
fit_RW_plus_WN <- function(model, pars.init = rep(log(var(diff(model$y))), 2)) {
objective <- function(pars) {
# Creiamo la funzione obiettivo con parametri ordinati come segue
# c(log(var(eta)), log(var(zeta)), log(var(epsilon)))
model$Q[1, 1, 1] <- exp(pars[1]) # prendendo l'exp le varianze saranno positive
model$H[1, 1, 1] <- exp(pars[2])
-logLik(model) # cambio segno perche optim() minimizza
}
# Minimizziamo la funzione obiettivo (usando l'algoritmo BFGS)
# e chiediamo di calcolare l'hessiana nel punto di massimo
optout <- optim(pars.init, objective, method = "BFGS", hessian = TRUE)
param_estim <- exp(optout$par) # stima di massima verosimiglianza delle varianze
mI <- solve(optout$hessian) # inversa della stima dell'info. di Fisher
mV <- diag(param_estim) %*% mI %*% diag(param_estim) # metodo delta in azione
sterr <- sqrt(diag(mV)) # vettore di errori standard
# Diamo i nomi ai parametri
names(param_estim) <- names(sterr) <-colnames(mV) <- rownames(mV) <- c("Q", "H")
list(hatparameters = param_estim, stderr = sterr, cov = mV, logLik = -optout$value, optim.out = optout)
}Vediamo ora la nostra funzione in azione.
fit1 <- fit_RW_plus_WN(modello1)
# valutiamo innanzitutto la convergenza: se 0 va bene
fit1$optim.out$convergence## [1] 0# guardiamo la stima dei parametri:
print(fit1$hatparameters)## Q H
## 1469.163 15098.651# stderr delle stime:
print(fit1$stderr)## Q H
## 1280.358 3145.560KFAS offre un metodo di stima dei parametri tramite la funzione fitSSM.
fitSSM(model, inits, updatefn) necessita di tre parametri formali:
La funzione fitSSM() restituisce una lista con optim.out, l’output della funzione optim() usata internamente per massimizzare la log-verosimiglianza, e model contenente il modello in forma state space con i parametri stimati (come abbiamo fatto nella funzione precedente).
updt <- function(pars, model) {
model$Q[1, 1, 1] <- exp(pars[1])
model$H[1, 1, 1] <- exp(pars[2])
model
}
# queste due sono uguali
# updt <- function(pars, model) {
# model$Q[[1]] <- exp(pars[1])
# model$H[[1]] <- exp(pars[2])
# model
# }
fit2 <- fitSSM(modello1, rep(log(var(diff(Nile))), 2), updt)
# valutiamo la convergenza dell'ottimizzazione: se 0 va bene
fit2$optim.out$convergence## [1] 0# Vediamo le stime nelle matrici
message("Q");fit2$model$Q[[1]];message("H");fit2$model$H[[1]]## Q## [1] 1466.32## H## [1] 15098.18# stderr delle stime:
fit3 <- fitSSM(modello1, rep(log(var(diff(Nile))), 2), hessian = TRUE)
V <- solve(fit3$optim.out$hessian)
GVG <- diag(exp(fit3$optim.out$par)) %*% V %*% diag(exp(fit3$optim.out$par))
stderr <- sqrt(diag(GVG))
message("Q stderr");stderr[1];message("H stderr");stderr[2]## Q stderr## [1] 1279.011## H stderr## [1] 3145.357Si noti che quando il parametro formale updatefn non viene specificato, la funzione fitSSM() usa come default una funzione che assume che i parametri da stimare siano le varianze lasciate a NA sulle diagonali delle matrici Q e H. Per garantire la positività delle varianze, tale funzione assume che i parametri in pars debbano essere trasformati con \(\exp(\theta)\).
fit3 <- fitSSM(modello1, rep(log(var(diff(Nile))), 2))
# valutiamo la convergenza dell'ottimizzazione: se 0 va bene
fit3$optim.out$convergence## [1] 0# Vediamo le stime nelle matrici
message("Q");fit3$model$Q;message("H");fit3$model$H## Q## , , 1
##
## [,1]
## [1,] 1466.32## H## , , 1
##
## [,1]
## [1,] 15098.18# stderr delle stime:
fit3 <- fitSSM(modello1, rep(log(var(diff(Nile))), 2), hessian = TRUE)
V <- solve(fit3$optim.out$hessian)
GVG <- diag(exp(fit3$optim.out$par)) %*% V %*% diag(exp(fit3$optim.out$par))
stderr <- sqrt(diag(GVG))
message("Q stderr");stderr[1];message("H stderr");stderr[2]## Q stderr## [1] 1279.011## H stderr## [1] 3145.357