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: