Numeric errors using Estimate with outputs with missing values

[Pedro Gea]

Parece que algo ha cambiado últimamente en torno al método Estimate pues la estimación de outputs con omitidos devuelve errores del tipo:

ERROR: [] No se pudieron calcular los valores de las interrupciones debido a un error de tipo numérico resolviendo el problema de programación lineal correspondiente mediante el método de descomposición de valor singular.
  Max(|PIW'(PIW*CI+A_)|)=0.212837490276996

que antes no devolvía.
Además, supongo que consecuencia de estos errores, los parámetros estimados ahora difieren de los que estimaba cuando no había errores.

Para comprobarlo se puede usar el código siguiente:

Serie omitido = Log((1-Pulse(y2000m03, C))*Exp(1));
Serie output = SubSer(Gaussian(10,1,C), y2000, y2001)* omitido;

Set Estimate(ModelDef(
  Serie output, 1, 0, 1, 0,  
  Polyn 1, SetOfPolyn(1-0.1*B), SetOfPolyn(1),
  Copy(Empty), Copy(Empty)
), y2000, y2001);

[Víctor de Buen Remiro]

Lo único que cambió, y de eso hace ya 6 meses en el ticket #774, es que se muestra un mensaje de error que antes no se mostraba, pero el error sí ocurría igualmente pero no se detectaba.

Lo que está ocurriendo es lo que dice el mensaje, que no se pueden calcular los omitidos, así de simple. Puede ser una carencia del método de cálculo pero no hay forma de introducir otro sin perder la eficiencia, aunque en todos los casos que he visto lo que ocurre es que el modelo no es viable porque hay un número excesivo de omitidos o bien tienen una estructura incompatible con el modelo ARIMA.

[Víctor de Buen Remiro]

Por otra parte, el que en una iteración intermedia se tropiece con ese problema no implica necesariamente que vaya a fracasar, puede ser que el estimador haya caído eventualmente en una zona de parámetros sin solución y luego se mueva a otra factible y acabe dando un modelo correcto.

Quizás se podría cambiar la categoría del mensaje de error a warning.

[Pedro Gea]

De acuerdo.

Entiendo que puede tratarse de alguna carencia del método de cálculo
ya que en algunos ejemplos que pruebo no parece haber nada inadecuado.
En ese caso, creo que debería ser una advertencia y no un error.

De aquí me surgen algunas dudas sobre cómo se estiman los valores omitidos
en el output, no sé si hay algo documentado sobre ello.

Por otro lado en el libro de Daniel Peña "Estadística. Modelos y métodos"
en la sección de "Regresión dinámica. 8. Análisis de intervención" encuentro
lo siguiente:
"Otra de sus utilizaciones [del análisis de intervención] es estimar valores
ausentes: si falta algún valor en la serie podemos sustituirlo por la media de
los dos adyacentes y estimar un modelo con una intervención en dicho punto.
Puede demostrarse (Peña, 1986) que si llamamos z_m al valor utilizado
para llenar el hueco y û a la estimación de la intervención, la predicción
optima del valor ausente es z_m - û y es independiente de z_m.
Este procedimiento puede usarse para completar la serie donde faltan algunas
observaciones introduciendo simultáneamente varias intervenciones."

Esto parece sugerir que los omitidos del output podrían estimarse poniendo
ceros donde aparezcan valores omitidos y añadiendo pulsos al conjunto de inputs.
¿Es esto así? ¿Se hacen en Estimate otras consideraciones?.

[Víctor de Buen Remiro]

Para ver si hay algún problema en el algoritmo necesitaría contar con un ejemplo convenientemente aislado de esos en los no parece haber nada inadecuado.

La documentación sobre el método la puedes ver en el capítulo 6.1 Interrupciones con continuidad inercial, del documento ​https://www.tol-project.org/export/HEAD/tolp/trunk/doc/modeling/ARIMA_Estimate/EstConInt2.doc EstConInt2.doc, y se corresponde al caso en el que el proceso realmente se interrumpe, no está definido para una fecha o intervalo temporal. En este caso el rellenado de la serie es completamente ficticio y se hace por motivos prácticos para poder aplicar los polinomios de la estructura ARIMA con lo que los valores introducidos deben ser los que menor discordancia introduzcan en las ecuaciones en diferencias.

En el caso del que habla Daniel Peña, la serie no se interrumpe en realidad, sino que por algún motivo un dato aislado, o unos pocos consecutivos, no son conocidos. Sólo en ese contexto tiene sentido aplicar la intervención con pulsos. La justificación se encuentra en el punto 6.2. En tal caso, el usuario de TOL debe filtrar los omitidos y añadir los pulsos correspondientes, lo cual le permite calcular los valores desconocidos. Por este motivo no está implementado internamente, pues se puede resolver redefiniendo el modelo y Estimate se reserva el concepto de interrupción para los valores omitidos que le lleguen. Obsérvese que en una misma serie output pueden aparecer omitidos de ambas especies y sólo el analista está en condiciones de distinguirlos, interviniendo unos y dejando otros.

Cuando los datos omitidos, bien por interrupción bien por desconocimiento conforman un intervalo demasiado largo, el método de cálculo debe ser otro (6.3), pues es necesario recalcular los valores iniciales de las ecuaciones en diferencias. Como eso es bastante complicado de implementar y además resultaría muy lento de ejecutar, lo que se hace en estos casos es modificar el fechado de la serie para extraerle el máximo un número de ciclos periódicos consecutivos para reducir el impacto de la discontinuidad inercial.

[Víctor de Buen Remiro]

(In [2118]) Refs #899

[Pedro Gea]

Está claro que la apreciación sobre los distintos tipos de omitidos, o al menos de las distintas maneras de estimarlos es muy interesante.

Aún no he tenido tiempo de leerlo con detenimiento, sin embargo haciendo algunas pruebas encuentro cosas inesperadas:

Adjunto el código de un "banco de pruebas", en el que es fácil variar el polinomio ARIMA y el número de omitidos seguidos de un modelo y estimarlo con las interrupciones de Estimate y con pulsos.
Observo dos cosas que me sorprenden:

• En el caso de un sólo omitido, los resultados con uno y otro método coinciden ( hasta el sexto decimal) con la antigua versión que tengo de TOLBase 1.1.7 (b.21 2009-10-09). Sin embargo con la versión actual, además de los errores que originaron este tique, esto no ocurre.
• En el caso de varios omitidos seguidos, el valor estimado para el omitido con el método de las interrupciones de Estimate es el mismo para todos, y claro, en esas circunstancias es lógico que no coincidan. ¿Es lógico que el valor estimado para varias interrupciones seguidas sea el mismo?

[Víctor de Buen Remiro]

Desde luego que no es lógico que salga el mismo valor para todos, ni tampoco que cambie el resultado de la versión pues en teoría sólo se añadió un mensaje de error.
Probaré el código para depurarlo.

Lo que sí estoy viendo es una cosa que está mal en el generador de ruido ARMA, aunque en este caso no tiene mucha importancia, pero que puede serlo en otras circunstancias. Veo que para inicializar la ecuación en diferencias del ruido usas normales independientes pero esto es incompatible con la estructura ARMA pues está claro que no pueden ser independientes si no se trata del modelo trivial o de ruido blanco, de grados AR(0)MA(0). Generar los valores iniciales tanto del ruido como de los residuos no es nada sencillo pero se puede aproximar bastante bien con el sencillo truco de alargar hacia atrás la serie de residuos de forma muy prolongada y efectuar la ecuación en diferencias con valores iniciales nulos tanto residuos como en ruido, y luego quedarse sólo con la parte del ruido que queda en el intervalo de interés. Si la prolongación hacia atrás fuera infinita el método sería exacto, aunque tardaría demasiado claro. El grado de aproximación para una prolongación finita dependerá del decaimiento exponencial de la estructura ARMA, por lo que con raíces cercanas a la unidad puede ser preciso un largo periodo. En general bastaría con tomar 100 veces el máximo de los grados AR y MA.

En el caso del ejemplo que propones sería algo así

Date begin = y2010;
Date end = y2010m12d31;

Serie residuals_gen = Gaussian(0, 1, C);

Polyn ar = 1 - 0.98*B; // - 0.2*B**2 +0.3*B**3;
Polyn ma = 1; // - 0.8*B;
Polyn dif = 1; // - B;

Real p = Degree(ar);
Real q = Degree(ma);
Real d = Degree(dif);
Real initLength = 100*Max(p+d,q);
Date begin_ = Succ(begin,Dating(residuals_gen),-initLength);
Serie residuals_= SubSer(residuals_gen,begin_,end);
Serie noise_ = DifEq(ma/(dif*ar), residuals_);
Serie residuals = SubSer(residuals_,begin, end);
Serie noise = SubSer(noise_, begin, end);

[Víctor de Buen Remiro]

(In [2127]) refs #899

[Pedro Gea]

Parece que es claro que Estimate no trata correctamente los omitidos, dejando a un lado las dudas sobre su significado: si son "interrupciones verdaderas" o "datos desconocidos".

En MMS permitimos la estimación de omitidos mediante Estimate.
No sé si el problema se podrá solucionar pronto o quizá deberíamos sustituir los omitidos por ceros y añadir pulsos para su estimación.

[Víctor de Buen Remiro]

Yo no sé lo que puede estar pasando porque incluso copiando directamente el código de CleanInterruptions del fichero modcalc.cpp de la versión 1.1.7 sigue fallando en la 2.0.1

Tiene que ser por tanto un efecto secundario de un cambio en otra parte del código al que se llame desde ahí.

Tendré que entrar a trazar en ambas versiones y ver dónde cambian los resultados. Se trata de un proceso de depuración tremendamente complicado y aunque estoy casi seguro de que lo encontraré no soy capaz de estimar el tiempo que puede llevar.

[Víctor de Buen Remiro]

El problema está en MinimumResidualsSolve que no sé qué le pasa pero falla en determinadas circunstancias. Como no veo una forma rápida de arreglarlo voy a a hacer un método basado en la descomposición de Cholesky de BVMat pues muy a menudo se trata de matrices sparse.
En cuanto haga las pruebas subiré los cambios

[Víctor de Buen Remiro]

(In [2193]) Fixes #899

[Víctor de Buen Remiro]

(In [2195]) Refs #899

[Víctor de Buen Remiro]

(In [2196]) Refs #899