Evaluación de pronóstico en R: validación cruzada de modelos de series temporales

En el siguiente ejemplo se muestra una validación cruzada de series temporales en R. Necesita evaluar la capacidad de pronóstico de un modelo para elegir alguna alternativa en función de su rendimiento. Es decir, necesito verificar cual hubiese sido el pronóstico histórico de una serie temporal, por ejemplo con un modelo AR(1), para comparar con la precisión de otros modelos y decidir cual estrategia seguir en lo adelante. En este sentido, necesitaríamos realizar el ejercicio de proyección en distintos puntos temporales, para posteriormente recuperar los errores y comparar estas medidas de error entre distintos modelos.

 1.       Simulamos la serie temporal (en el caso de una base de datos, deberías importarlas en R).

 library(writexl)   # Exportar a Excel
library(readxl)
library(lubridate) # dates
library(tidyverse)
library(forecast)
library(timeSeries)
library(zoo)
 
# Importando datos ---------------
# ********************************************************************
# ********************************************************************
 
# Configuración inicial
set.seed(4)    # Asegura reproducibilidad
n <- 250       # Número de días (simulación de aproximadamente un año laboral)
mu <- 0        # Promedio del cambio diario (sin tendencia)
sigma <- 0.01  # Volatilidad diaria
 
# Generación de la caminata aleatoria
tipo_cambio <- cumsum(c(1, rnorm(n, mean = mu, sd = sigma)))  # Inicia en 100 y suma los cambios diarios
 
# Creación de un data.frame con fechas
fechas <- seq.Date(from = as.Date("2023-01-01"), by = "day", length.out = n + 1)
df <- data.frame(fecha = fechas, tipo_cambio = tipo_cambio)
 
# Gráfico de series temporales
library(ggplot2)
ggplot(df, aes(x = fecha, y = tipo_cambio)) +
  geom_line(color = "blue") +
  labs(title = "Simulación de Caminata Aleatoria del Tipo de Cambio",
       x = "Fecha",
       y = "Tipo de Cambio") +
  theme_minimal()
 
library(xts)
tipo_cambio_xts <- xts(tipo_cambio, order.by = fechas)

 
2.       Validación cruzada

En resumen, lo que se hace es definir una ventana de estimación y otra de validación (pronóstico). En este ejemplo son de 60 y 5 observaciones respectivamente. 60 observaciones serán usadas para estimar el modelo, mientras que 5 será la longitud de la ventana del pronóstico. De esta forma, en la primera iteración se va estimar los modelos con observaciones del 1:60 y otra ventana de pronóstico de la 61-65, luego 2:61 y 62:66, … y así respectivamente hasta recorrer toda la muestra. El documento de trabajo “Pronóstico del consumo privado en la República Dominicana” de Alejandro J. Balcácer y Nerys F. Ramírez, muestra en detalle esta estrategia.

 # Análisis de precisión. Error histórico ---------------

h.valid <- 5
horizonte <- paste0("h", 1:h.valid)
vent <- 60

Posteriormente, se genera un for para realizar este recorrido usando una serie. Primero se define el valor del ultimo dato observado en cada iteración del bucle (seq_for), este dato sirve de referencia para el resto de las ventanas. Además, se define un contador para el punto de inicio de la ventana de estimación (cont). Los puntos clave de esta secuencia son:

-          Indicar la secuencia con la posición de las observaciones que vamos a usar en cada punto de la serie temporal [obser in seq_for].

-          La serie temporal se segmenta en dos partes: la estimación (tipo_cambio_xts[cont:obser]) y la otra parte de pronósticos (tipo_cambio_xts[(obser+1):(obser+h.valid)]). Fíjese que la primera va desde el contador hasta obser, recuerde que obser es el valor final hasta donde irá la ventana de estimación y el valor que en el presente ejemplo sirve de referencia. Luego, la parte para evaluar las proyecciones se definen a partir de la última observación de la parte de estimación (obser+1):(obser+h.valid).

-          La tercera parte importante del bucle es incluir el modelo de pronóstico. Aquí debemos usar alguna estrategia para proyectar los periodos que nos resulten de interés. En este ejemplo son 5 observaciones, y se usa un modelo ar(1).

-          Finalmente, se crea un data.frame para guardar el pronóstico, el dato observado, los horizontes y el nombre del model. Además, se guarda cada data.frame en una lista creada anteriormente (data_Rbind). Primero se llaman todos los objetos que inician con “df_” (tener pendiente guardar cada nuevo modelo con ese nombre).  

 # Esto lo que hace, es que permite recorrer las ventanas
seq_for <- vent:(nrow(tipo_cambio_xts)-h.valid)
 
# OJO CORRER EL BULCE DESDE AQUI
cont <- 1
data_Rbind <- list() # Lista que guarda resultados
 
for (obser in seq_for){
  
  tc_train <- tipo_cambio_xts[(obser-vent+1):obser]
  tc_test <- tipo_cambio_xts[(obser+1):(obser+h.valid)]
 
  fecha_ciclo <- fechas[obser]
  fecha_hat <-  fechas[(obser+1):(obser+h.valid)]
 
ordena_data <- function(objetos,...){
  # objetos = lista de nombre de vectore
 
  # Iterar sobre los nombres de los objetos y cambiar los nombres de las columnas
  for (i in seq_along(objetos_df)) {
    xdf_name <- get(objetos_df[i])
    new_column_names <- c("fechaCiclo","fechaHat","horizonte","obs","pro","model")
    names(xdf_name) <- new_column_names
    assign(objetos_df[i], xdf_name)
  }
 
  lista_objetos <- lapply(objetos_df, get)
  df1_final <- dplyr::bind_rows(lista_objetos)
  return(df1_final)
}
    ### 1. ESTRATEGIA AR(1) . Ventana móvil ************************
 
    model_ar1 <- Arima(tipo_cambio_xts, order=c(1,0,0))
    ht_f_ar1 <- forecast::forecast(model_ar1, h=h.valid)$mean                 
   
    df_ar1 <- data.frame(fecha_ciclo, fecha_hat, horizonte,  #fechas
                         obsT=tc_test, # dato observado
                         proT=ht_f_ar1,
                         model="1.Ar(1)_móvil")
 
    ### 1. ESTRATEGIA AR(2,1) . Ventana móvil ************************
   
    model_ar21 <- Arima(tipo_cambio_xts, order=c(2,0,1))
    ht_f_ar21 <- forecast::forecast(model_ar21, h=h.valid)$mean                 
   
    df_ar21 <- data.frame(fecha_ciclo, fecha_hat, horizonte,  #fechas
                         obsT=tc_test, # dato observado
                         proT=ht_f_ar21,
                         model="2.Ar(21)_móvil")
   
### ******************************************************************
### ******************************************************************
objetos <- ls()
objetos_df <- objetos[grep("^df_", objetos)]
 
data_Rbind[[(obser-vent+1)]] <- ordena_data(objetos_df)
   
# Esto al final siempre       
cont <- cont+1  # cont y esto es equivalente (obser-vent+1):
}

Luego de culminar, se utiliza una función propia para ordenar la data (ordena_data) y convertirla en un único data.frame con toda la data:

data_forecast_final
# A tibble: 1,870 × 6
   fechaCiclo fechaHat   horizonte   obs  proy model        
   <date>     <date>     <chr>     <dbl> <dbl> <chr>        
 1 2023-03-01 2023-03-02 h1        0.969 1.01  1.Ar(1)_móvil
 2 2023-03-01 2023-03-03 h2        0.981 1.01  1.Ar(1)_móvil
 3 2023-03-01 2023-03-04 h3        0.972 1.01  1.Ar(1)_móvil
 4 2023-03-01 2023-03-05 h4        0.978 1.00  1.Ar(1)_móvil
 5 2023-03-01 2023-03-06 h5        0.987 1.00  1.Ar(1)_móvil
 6 2023-03-01 2023-03-02 h1        0.969 1.01  1.Ar(21)_móvil
 7 2023-03-01 2023-03-03 h2        0.981 1.01  1.Ar(21)_móvil
 8 2023-03-01 2023-03-04 h3        0.972 1.01  1.Ar(21)_móvil
 9 2023-03-01 2023-03-05 h4        0.978 1.00  1.Ar(21)_móvil
10 2023-03-01 2023-03-06 h5        0.987 0.999 1.Ar(21)_móvil

 Con esta data, ahora se puede realizar cualquier análisis sobre los errores promedio, su evolución histórica o el comportamiento de estos errores condicionados a eventos económicos, por ejemplo, los errores durante el COVID o periodos de alta incertidumbre.

 # Cuadrando la base de los errores
# Crea la data final con todos los los dtos necesarios
data_forecast_final <- dplyr::bind_rows(data_Rbind) |> as_tibble()
names(data_forecast_final) <- c("fechaCiclo","fechaHat","horizonte","obs","proy","model")
 
# Errores total
# Total
data_forecast_final |>
  #dplyr:: filter(fechaCiclo > as.Date("2018-01-01"))  |>
  mutate(error=obs-proy) |>
  group_by(horizonte, model) |> # (horizonte, model, COVID)
  summarize(
    mae_e_abs = mean(abs(error), na.rm = TRUE),
  ) |>
  pivot_wider(names_from = horizonte, values_from = mae_e_abs) |>
  arrange(h1) |>
  dplyr::select(model,horizonte, everything())
 
`summarise()` has grouped output by 'horizonte'. You can override using the `.groups` argument.
# A tibble: 2 × 6
  model              h1     h2     h3     h4     h5
  <chr>           <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 1.Ar(21)_móvil 0.0479 0.0442 0.0405 0.0377 0.0350
2 1.Ar(1)_móvil  0.0483 0.0443 0.0409 0.0378 0.0353

Modelos de frecuencia mixta en R

 1.    Simulando la data

La siguiente entrada muestra un ejemplo de como estimar un modelo de frecuencia mixta. Puntualmente simulamos una serie de consumo trimestral que queremos estimar en una regresión contra una serie de ventas mensuales. Como tal, la diferencia en el numero de observaciones limita la posibilidad de estimar directamente un modelo de regresión tradicional -aquí suponemos que las series son estacionarias, proceso obligatorio para estimar una regresión entre series temporales-. Todo el ejercicio se realiza en Package ‘midasr’ (Zemlys, 2022).

En primer lugar, simulamos la serie de consumo:

# Simulando series
library(ggplot2)
library(forecast)
library(midasr)
 
# Simulando consumo
# Definición de parámetros
set.seed(1234)
trend <- 0.5  # Tendencia de la serie
seasonality <- 0.2  # Estacionalidad trimestral
n <- 22*4  # Número de trimestres a simular (22 años)
 
# Generación de la serie trimestral
consumo <- rnorm(n, mean = trend + seasonality * rep(1:4, length.out = n))
consumo_trimestral <- ts(consumo, frequency = 4, start = 2000)
 
autoplot(consumo_trimestral)

 

Posteriormente simulamos la serie de ventas, teniendo presente que el número de observaciones de la variable mensual, debe coincidir con las necesarias para equiparar las observaciones de mi serie en frecuencia más alta. Puntualmente, si tenemos una serie trimestral de 10 años, entonces necesitaríamos 4*10 observaciones de mi serie en menor frecuencia.

# simulando ventas
# Definición de parámetros
trend <- 0.8  # Tendencia de la serie
seasonality <- c(0.5, 0.3, 0.2)  # Estacionalidad mensual (ejemplo: enero, febrero, marzo)
n2 <- (n/4)*12  # Número de meses a simular
 
# Generación de la serie mensual
ventas_mensuales <- numeric(n)
for (i in 1:n2) {
  ventas_mensuales[i] <- rnorm(1, mean = trend + seasonality[(i - 1) %% length(seasonality) + 1])
}
 
ventas_mensuales <- ts(ventas_mensuales, frequency = 12, start = 2000)
 
autoplot(ventas_mensuales)

2.     Estimando el modelo midas

Ahora que hemos simulado nuestras series, podemos estimar nuestro modelo en frecuencia mixta. Primero utilizamos una ventana de datos para poder comparar los valores proyectados con los valores observados. Luego usamos la función midas_r para estimar el modelo en frecuencia mixta. Note aquí que la clase está en la función fmls que realiza un collapse de la data en alta frecuencia para que este coincida en el total de observaciones. Dentro de esta función se coloca la variable independiente, el numero de rezagos a incluir, la cantidad de veces que esta contenida mi serie de baja frecuencia en la de alta (trimestral y mensual 3; anual y mensual 12…)…

 # train data
x_q <-  window(ventas_mensuales, end = c(2019, 12))
y_c <-  window(consumo_trimestral, end = c(2019, 4))
aa <- round(length(x_q)/length(y_c))
 
trend <- 1:length(y_c)
 
mr <- midas_r(y_c ~ trend + fmls(x_q, 11, aa, nealmon), start = list(x_q = rep(0, 3)))
 
MIDAS regression model with "ts" data:
Start = 2000(4), End = 2019(4)
 model: y_c ~ trend + fmls(x_q, 11, aa, nealmon)
(Intercept)       trend        x_q1        x_q2        x_q3
    0.50057     0.01001    -0.12472    20.67112    -0.84271
 
Function optim was used for fitting

 [Nota: note que la funciones fmls selecciona volares puntuales de una vector de datos en el caso de colocar 0 toma la observación actual. 1, toma la actual y crea un vector del tamaño de la serie de baja frecuencia que adhiere el vector anterior. 2, a los dos vectores anteriores, agrega uno con el segundo rezago. Finalmente, el último argumento le indica cada cuantas observaciones de tu serie de alta frecuencia tu tomaras un dato para crear un vector cuya longitud coincida con el de baja frecuencia.

x<-c(1:12)
x
 [1]  1  2  3  4  5  6  7  8  9 10 11 12
 
fmls(x, 0, 3)
     X.0/m
[1,]     3
[2,]     6
[3,]     9
[4,]    12
 
fmls(x, 1, 3)
     X.0/m X.1/m
[1,]     3     2
[2,]     6     5
[3,]     9     8
[4,]    12    11
]

Ahora, solo témenos que considerar cuales son los valores futuros de las variables independientes. Esto con el fin de poder crear pronósticos condicionados.

##New trend values
n.valid <- 8
xn <- forecast::snaive(x_q,aa*n.valid)$mean
trendn <- length(y_c) + 1:n.valid

Finalmente, usamos la función forecast para obtener los valores predichos para mi variable independiente.

fh_midas_1 <- forecast(mr,  list(trend = trendn, x_q = xn), method = "static")
fh_midas_1
         Qtr1     Qtr2     Qtr3     Qtr4
2020 1.213594 1.255759 1.216617 1.158818
2021 1.253628 1.295793 1.256651 1.198851

 3.     Estimando un modelo dinámico (inluyendo un ar de y)

 En el siguiente ejemplo usamos mls para ingresar en el modelo valores rezagados de la variable dependiente.

 mr.dyn <- midas_r(y_c ~ trend +
                    mls(y_c, 1:2, 1, "*") +
                    fmls(x_q, 12, aa, nealmon),
                    start = list(x_q = rep(0, 3)))
 
MIDAS regression model with "ts" data:
Start = 2001(3), End = 2019(4)
 model: y_c ~ trend + mls(y_c, 1:2, 1, "*") + fmls(x_q, 18, aa, nealmon)
(Intercept)       trend        y_c1        y_c2        x_q1        x_q2        x_q3
   0.531114    0.009247    0.105423   -0.017154   -0.214118    2.527028   -0.111902
 
Function optim was used for fitting
 
fh_midas_2 <- forecast(mr.dyn,  list(trend = trendn, x_q = xn), method = "dynamic")
fh_midas_2
         Qtr1     Qtr2     Qtr3     Qtr4
2020 1.223314 1.274915 1.186098 1.196951
2021 1.278432 1.315320 1.226402 1.237496

 4.     Comparando diversas estrategias de pronósticos

 Finalmente armamos una base de datos, un poco tosca porque seguro hay formas de hacer esto de manera mas eficiente y elegante. Adicionalmente, agregamos una alternativa de estimación mediante el modelo auto.arima. es importante verificar que mayor complejidad no siempre significa mejor precisión en la estimación de pronósticos.

 # AutoArima
fh_arima <- forecast::forecast(auto.arima(y_c, stepwise = FALSE))
 
### Graficando
obs_data <- window(consumo_trimestral, start = c(2020, 1))
 
data_fores <- data.frame(
  y_obs = c(y_c, rep(NA, n.valid)),
  f.midas1 = c(rep(NA, length(y_c)-1), tail(y_c,1), fh_midas_1$mean),
  f.midas2 = c(rep(NA, length(y_c)-1), tail(y_c,1), fh_midas_2$mean),
  fh_arima = c(rep(NA, length(y_c)-1), tail(y_c,1), fh_arima$mean),
  obs = c(rep(NA, length(y_c)-1), tail(y_c,1), obs_data),
  fecha = time(consumo_trimestral)
)
 
data_fores |> tail(8*4) |>
ggplot() +
  geom_line(aes(fecha,y_obs,colour="Obs"), size=1) +
  geom_line(aes(fecha, obs,colour="Obs1"),linetype = 2, size=1) +
  geom_line(aes(fecha, f.midas1,colour="xMidas"),linetype = 3, size=1) +
  geom_line(aes(fecha, f.midas2,colour="xArMidas"), linetype = 3, size=1) +
  geom_line(aes(fecha, fh_arima,colour="Arima"),linetype = 3, size=1) +
  theme_classic() +
  scale_color_manual(name = "Serie", values = c("Obs" = "black", "Obs1" = "black", "xMidas" = "darkred", "xArMidas" = "darkblue", "Arima"="darkgreen")) +
  theme(legend.position = "bottom")

Referencias

Zemlys, V. (2022). Mixed Data Sampling Regression. Package ‘midasr’. R CRAN.

Pronóstico directo vs. iterado: metodología y validación cruzada

En la siguiente entrada se desarrolla una función de R para proyectar paso a paso en un horizonte determinado. Puntualmente, se comparan dos estrategias: una donde se usa la función auto.arima (ver) para proyectar todos los valores en un horizonte determinado a partir del modelo identificado; y una segunda estrategia que consiste en proyectar a un solo paso, tomar este valor proyectado como observado y posteriormente volver a identificar el modelo ARIMA correspondiente. Posteriormente, se repite este procedimiento tantas veces como sea especificado en la ventana de pronósticos.

 Para el ejemplo utilizaremos la serie temporal AirPassengers:

 library(dplyr)
library(readxl)
library(tibble)
library(ggplot2)
library(forecast)
library(tseries)
require(lubridate)
library(tidyr)
 
autoplot(AirPassengers) +
  theme_minimal()

Ahora definimos el horizonte de pronóstico y usamos la funciones auto.arima y forecast (como hay varios paquetes que utilizan esta función, se indexa dentro del paquete para hacer referencia a esta y evitar ambigüedades en otros programas) para obtener las proyecciones del modelo ARIMA seleccionado por el algoritmo de automatización:

h.valid <- 12
 
# alternativa 1: direct
for1 <- forecast::forecast(auto.arima(AirPassengers), h=h.valid)
for1
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Jan 1961       445.6349 430.8903 460.3795 423.0851 468.1847
Feb 1961       420.3950 403.0907 437.6993 393.9304 446.8596
Mar 1961       449.1983 429.7726 468.6240 419.4892 478.9074
Apr 1961       491.8399 471.0270 512.6529 460.0092 523.6707
May 1961       503.3945 481.5559 525.2330 469.9953 536.7937
Jun 1961       566.8624 544.2637 589.4612 532.3007 601.4242
Jul 1961       654.2602 631.0820 677.4383 618.8122 689.7081
Aug 1961       638.5975 614.9704 662.2246 602.4630 674.7320
Sep 1961       540.8837 516.9028 564.8647 504.2081 577.5594
Oct 1961       494.1266 469.8624 518.3909 457.0177 531.2356
Nov 1961       423.3327 398.8381 447.8273 385.8715 460.7939
Dec 1961       465.5076 440.8229 490.1923 427.7556 503.2596

Ahora utilizamos la función propia AutoArima_iterated para realizar el procedimiento descrito anteriormente. En este sentido se hace una proyección a un paso [forecast::forecast(for_x, h=1)], luego esta proyección se toma como observada [x<-c(x[-1],fhat)], aunque eliminamos la primera observación del vector x [x[-1]] esto para evitar la ventana se mueva de forma recursiva (creciendo una observación en cada paso), y se mueva entonces como una ventana móvil, manteniendo el tamaño de la ventana durante todo el horizonte. Una vez se modifica este valor observado, obtenemos un nuevo vector de datos observados para alimentar el modelo, pero teniendo el dato que proyectamos anteriormente como valor observado. Este proceso se repite n veces, como cuantos horizontes deseamos obtener.

# alternativa 2: iterated
source("AutoArima_iterated.R")
 
function(x,hn=h.valid) {
  for (hj in 1:hn){
    for_x <- auto.arima(x)
    fhat <- forecast::forecast(for_x, h=1)$mean
    x<-c(x[-1],fhat)  # -1 mantiene el tamano de la ventana
  }
  return(tail(x,h.valid))
}
 
for2 <-AutoArima_iterated(x=AirPassengers,hn=h.valid)
for2
 
       Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov
1961 445.6349 460.7414 470.6938 497.8221 502.7595 495.1462 484.4825 495.1744 496.0780 491.9802 487.1089
          Dec
1961 483.7122

Ahora, para representar ambas proyecciones en un gráfico usando la función autoplot, copiamos las propiedades guardada en la serie temporal del primer pronostico (for1) en este segundo vector de datos (for2). En este caso, recordemos que la función original de pronóstico es forecast y esta guarda el pronóstico en formato de serie temporal, mientras que el segundo lo hace como un vector.

 for2 <- ts(for2, star=start(for1$mean), frequency = frequency(for1$mean))
 
tail(AirPassengers,24) %>%
autoplot() +
  autolayer(for1$mean, color = "blue") +
  autolayer(for2, color = "red") +
  theme_minimal()

Finalmente, es relevante realizar una validación cruzada de estas proyecciones. Esto es principalmente útil en el marco de un sistema de pronósticos, pues nos ayuda a evaluar en el tiempo cual de las estrategias ha mostrado mayor precisión (menor error de pronósticos). En tal sentido usamos la función tsCV del paquete forecast de R dato que permite fácilmente obtener una serie histórica de los errores de pronóstico a cada horizonte.

1.       Primero definimos una función para proyectar valores (far1).

2.       Luego utilizamos la función de validación cruzada tsCV para obtener un histórico de errores a cada horizonte. Como son 6, aquí obtendremos una base de datos con 6 columnas y tantas observaciones como tenga nuestra data original.

3.       Finalmente obtenemos el RMSE de estos errores. Raíz del promedio de errores al cuadrado.

Sqrt(mean(e^2)) 

far1<- function(x, h){forecast::forecast(auto.arima(x), h=h)}
eAr1 <- tsCV(AirPassengers, forecastfunction = far1, h = h.valid, window=vent)
e_Ar1 <- sqrt(apply(eAr1^2, MARGIN=2, FUN=mean, na.rm = TRUE))

El procedimiento anterior lo repetimos, pero ahora usando la función AutoArima_iterated. Note que la única diferencia respecto a la evaluación anterior, es que usamos la función tsCv2, esta es la misma función tsCV, pero modificada para aceptar trabajar con cualquier pronóstico, dado que la función original solo trabaja con pronósticos arrojado en la estructura de los pronósticos del paquete forecast. Esto constituye una limitante para evaluar otras estrategias cuyo formato no coincidan. Ahora se modifica la función para que trabaje con cualquier pronostico siempre que coincida la longitud de los mismos. Dicha función modificada la colocamos como anexo en este documento.

# Auto.Arima paso a paso
source("tsCv2.R")
 
eAr1b <- tsCv2(AirPassengers, forecastfunction = AutoArima_iterated, h = h.valid, window=vent)
e_Ar1b <- sqrt(apply(eAr1b^2, MARGIN=2, FUN=mean, na.rm = TRUE))
e_Ar1b

 Ahora comparamos ambas estrategias, verificándose que la segunda estrategia es menos adecuada al intentar anticipar la serie, especialmente porque no parece capturar muy bien el componente estacional de la serie:

 t(rbind(e_Ar1,e_Ar1b))
    e_Ar1   e_Ar1b
h=1  13.42502 13.42502
h=2  16.01273 43.62504
h=3  18.80455 68.52383
h=4  20.42555 81.19967
h=5  20.51924 86.99055
h=6  20.98445 85.87618
h=7  21.54783 82.59922
h=8  22.37709 81.59745
h=9  22.52545 77.86818
h=10 23.80513 72.09909
h=11 24.62540 62.98470
h=12 25.83688 48.98253

Anexo:

tsCv2 <- function (y, forecastfunction, h = 1, window = NULL, xreg = NULL,
          initial = 0, ...)
{
  y <- as.ts(y)
  n <- length(y)
  e <- ts(matrix(NA_real_, nrow = n, ncol = h))
  if (initial >= n)
    stop("initial period too long")
  tsp(e) <- tsp(y)
  if (!is.null(xreg)) {
    xreg <- ts(as.matrix(xreg))
    if (NROW(xreg) != length(y))
      stop("xreg must be of the same size as y")
    xreg <- ts(rbind(xreg, matrix(NA, nrow = h, ncol = NCOL(xreg))),
               start = start(y), frequency = frequency(y))
  }
  if (is.null(window))
    indx <- seq(1 + initial, n - 1L)
  else indx <- seq(window + initial, n - 1L, by = 1L)
  for (i in indx) {
    y_subset <- subset(y, start = ifelse(is.null(window),
                                         1L, ifelse(i - window >= 0L, i - window + 1L, stop("small window"))),
                       end = i)
   
      fc <- forecastfunction(y_subset, h = h)
   
      e[i, ] <- y[i + seq(h)] - fc[seq(h)]
  }
  if (h == 1) {
    return(e[, 1L])
  }
  else {
    colnames(e) <- paste("h=", 1:h, sep = "")
    return(e)
  }
}