pronosticoNP
Pronóstico No Paramétrico en R
Rodríguez y Siado (2003) ilustran una forma no paramétrica de realizar pronósticos de series de tiempo estacionarias, empleando un estimador kernel gaussiano para la media condicional. A continuación se encuentra un procedimiento en R para realizar este tipo de pronósticos.
Estimador kernel:
rn = function(z,m=1,d=1,hn=.5){
n = length(z) ; Xn = z[n:(n-d+1)] ; Y = z[(d+m):n]
X = numeric(0) ; for(t in d:(n-m)) X = rbind(X,z[t:(t-d+1)])
Kdt = apply(X,1,function(x) (2*pi)^(-d/2)*exp(-sum(((Xn-x)/hn)^2/2)))
sum(Kdt*Y/sum(Kdt))
}
MELM (calidad del pronóstico):
MELM = function(hn,d,m,z){
n = length(z) ; p = ifelse(n<100,floor(n/4),floor(n/5)) ; melm = 0
for(r in 0:(p-m)) melm = melm+abs(rn(z[1:(n-p+r)],m=m,d=d,hn=hn)-z[(n-p+r+m)])
melm/(p-m+1)
}
Pronóstico:
sea:
m: horizonte de pronóstico. Por defecto m=1.
d: vector de coeficientes de Markov para cada m. Es un vector de longitud m.
hn: ancho de banda, 0<hn<1. Si hn=NULL, entoces para cada horizonte 1:m la función busca el hn que minimiza el MELM
dmax: si d=NULL, entonces la función hace d = 1:dmax. Por defecto es 10.
predNP = function(z,m=1,d=NULL,hn=NULL,dmax=10){
m = 1:m
if(!is.null(d)) dmh = data.frame(d,m)
if(is.null(d)){
d = 1:dmax ; grid = expand.grid(d=d,m=m)
op = apply(grid,1,function(x) unlist(optimize(MELM,c(0,1),d=x[1],m=x[2],z=z)))
dmh = data.frame(grid,t(op))
dmh = dmh[is.element(dmh[,4],tapply(dmh[,4],dmh[,2],function(x) min(x))),]
}
pred = apply(dmh,1,function(x) rn(z,m=x[2],d=x[1],hn=x[3]))
data.frame(dmh,pred)
}
Como ejemplo, realicemos un pronóstico de 12 meses para la serie AirPassengers de R:
x = log10(AirPassengers) # Serie transformada
z1 = diff(x,lag=12) # Diferencia estacional
z = diff(z1) # Diferencia no estacional
z.predNP = predNP(z,12) # Pronóstico no paramétrico de la serie transformada y diferenciada
z.completa = c(z,z.predNP$pred) # Serie completa
inv1 = diffinv(z.completa,lag=1,xi=z1[1]) # Inversa de la 1a diferencia
inv2 = diffinv(inv1,lag=12,xi=x[1:12]) # Inversa de la diferencia estacional
NPpred = ts(10^(inv2[-(1:144)]),freq=12,start=c(1961,1)) # Transformada inversa del pronóstico
NPpred
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1961 444.5 417.8 450.7 490.8 501.6 562.2 654.9 637.7 535.7 484.3 412.7 461.3
fit = arima(log10(AirPassengers),c(0,1,1),seasonal=list(order=c(0,1,1),period=12)) # Modelo SARIMA para comparar
pred = 10^predict(fit, n.ahead = 12)$pred # Pronóstico del modelo SARIMA
pred
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1961 450.4 425.7 479.0 492.4 509.1 583.3 670.0 667.1 558.2 497.2 429.9 477.2
ts.plot(AirPassengers, pred, NPpred, col = c(1,2,4), lwd=2) # Gráfico del pronóstico
legend(locator(1),c("SARIMA","No Paramétrico"),col=c(2,4),lty=1,lwd=2)# Click sobre el gráfico para agregar la leyenda
Comparte esto: