Introducción

La publicación de la obra Time Series Analysis: Forecasting and Control por G. E. P. Box y G. M. Jenkins en 1976 estableció un punto de inflexión en las técnicas cuantitativas de predicción en Economía. La metodología propuesta por estos autores, también conocida como metodología ARIMA, trata de realizar previsiones acerca de los valores futuros de una variable utilizando únicamente como información la contenida en los valores pasados de la propia serie temporal. Este enfoque supone una alternativa a la construcción de modelos uniecuacionales o de ecuaciones simultáneas, pues supone admitir que las series temporales poseen un carácter estocástico, lo que implica que deben analizarse sus propiedades probabilísticas para que éstas “hablen por sí mismas”.

El análisis univariante de series temporales presenta como ventaja frente a otros métodos de predicción el no depender de los problemas de información asociados a las variables endógenas o exógenas. Como hemos visto en capítulos anteriores, los modelos económicos que hemos estimado hasta el momento requerían un conjunto de variables exógenas que se utilizaban para explicar el comportamiento de una variable endógena. Sin embargo, en muchas ocasiones no se dispone de observaciones para alguna de las variables exógenas, ya sea porque no es posible medir la variable (por ejemplo, las expectativas de los agentes) o porque la muestra de datos de que disponemos para representar dicha variable presenta errores de medida. Este problema desaparece cuando se trata de modelizar una variable endógena mediante un modelo de tipo univariante como el propuesto por Box y Jenkins, ya que se hace depender a dicha variable tan sólo de su propio pasado y un conjunto de perturbaciones aleatorias, pero no de otras variables, caracterizando así las series económicas en su dimensión temporal.

En el presente apartado vamos a definir y caracterizar una amplia familia de estructuras estocásticas lineales, así como la metodología a seguir para seleccionar aquel modelo univariante que resulte más adecuado para representar la estructura estocástica de la variable económica que estemos analizando.

Procesos Estocásticos

Podemos definir un proceso estocástico como un conjunto de variables aleatorias asociadas a distintos instantes del tiempo. Así, en cada período o momento temporal se dispone de una variable que tendrá su correspondiente distribución de probabilidad; por ejemplo, si consideramos el proceso \(Y_t\), para \(t = 1\), tendremos una variable aleatoria, \(Y_1\), que tomará diferentes valores con diferentes probabilidades.

La relación existente, por tanto, entre una serie temporal y el proceso estocástico que la genera es análoga a la que existe entre una muestra y la población de la que procede, de tal forma que podemos considerar una serie temporal como una muestra o realización de un proceso estocástico, formada por una sola observación de cada una de las variables que componen el proceso. La tarea del investigador será, por tanto, inferir la forma del proceso estocástico a partir de las series temporales que genera.

Un proceso estocástico, \(Y_t\), se suele describir mediante las siguientes características: esperanza matemática, varianza, autocovarianzas y coeficientes de autocorrelación.

La esperanza matemática de \(Y_t\) se traduce en la sucesión de las esperanzas matemáticas de las variables que componen el proceso a lo largo del tiempo, tal que:

\[E(Y_t)=\mu_t,\ \ t=1,2,3...\]

Por su parte, la varianza de un proceso aleatorio es una sucesión de varianzas, una por cada variable del proceso:

\[Var(Y_t)=E(Y_t–\mu_t)^2,\ \ t=1,2,3...\]

Las autocovarianzas, por su parte, son las covarianzas entre cada par de variables del proceso, tales que:

\[\gamma_k=Cov(Y_t,Y_{t+k})=E[(Y_t–\mu_t)(Y_{t+k}–\mu_{t+k})]=\gamma_{t,t+k},\ \ t=1,2,3...\]

Finalmente, los coeficientes de autocorrelación son los coeficientes de correlación lineal entre cada par de variables que componen el proceso:

\[\rho_{t,t+k}=\frac{\gamma_{t,t+k}}{\sqrt{Var(Y_t)Var(Y_{t+k})}},\ \ t=1,2,3...,\ \ -1 \leq \rho_{t,t+k} \leq 1\]

Por último, a partir de los coeficientes de autocorrelación, vamos a definir dos funciones que nos serán muy útiles a lo largo del presente tema:

• Por un lado, la función de autocorrelación simple (fas) o correlograma, la cual es la representación gráfica de los coeficientes de autocorrelación en función de los distintos retardos o desfases entre las variables.

• La función de autocorrelación parcial (fap), que mide la correlación existente entre dos variables del proceso en distintos períodos de tiempo, pero una vez eliminados los efectos sobre las mismas de los períodos intermedios. Por ejemplo, puede que exista cierta correlación entre \(Y_t\) e \(Y_{t-2}\), debido a que ambas variables estén correlacionadas con \(Y_{t-1}\).

Dado que en la práctica se dispone de una muestra de un proceso estocástico, \(Y_1, ...Y_T\), se pueden obtener los coeficientes de autocorrelación simple y parcial muestral y utilizarlos como estimadores de los parámetros de la función de autocorrelación simple y parcial teórica.

Así, la función de autocorrelación simple (fas) puede estimarse a partir de las autocovarianzas del proceso tal que:

\[\hat \rho_{k}=\frac{\hat \gamma_k}{\hat \gamma_0}\]

Siendo:

\[\hat \gamma_0=\frac{\sum_{t=1}^{T}(Y_t-\bar Y)^2}{T}\]

\[\hat \gamma_k=\frac{\sum_{t=k+1}^{T}(Y_t-\bar Y)(Y_{t+k}-\bar Y)}{T-k}\]

La estimación de los parámetros de la función de autocorrelación parcial (fap) resulta algo más compleja, por lo que se verá en epígrafes posteriores.

Procesos estacionarios

Se dice que un proceso estocástico es estacionario en sentido estricto si todas las variables aleatorias que componen el proceso están idénticamente distribuidas, independientemente del momento del tiempo en que se estudie el proceso.

Es decir, la función de distribución de probabilidad de cualquier conjunto de k variables (siendo k un número finito) del proceso debe mantenerse estable (inalterable) al desplazar las variables s períodos de tiempo tal que, si \(P(Y_{t+1}, Y_{t+2}, …, Y_{t+k} )\) es la función de distribución acumulada de probabilidad, entonces:

\[P(Y_{t+1},Y_{t+2},…,Y_{t+k}) = P(Y_{t+1+s}, Y_{t+2+s},…,Y_{t+k+s}),\ \ \forall t, k, s\]

Sin embargo, la versión estricta de la estacionariedad de un proceso suele ser excesivamente restrictiva para las necesidades prácticas de un economista. Es por ello que generalmente nos conformaremos con un concepto menos exigente, el de estacionariedad en sentido débil o de segundo orden, la cual se da cuando la media del proceso es constante e independiente del tiempo, la varianza es finita y constante, y el valor de la covarianza entre dos periodos depende únicamente de la distancia o desfase entre ellos, sin importar el momento del tiempo en el cual se calculan. Dicho de otro modo, todos los momentos de primer y segundo orden de un proceso estocástico que sea estacionario en sentido débil deben ser invariantes en el tiempo.

La contrastación empírica de algunas de estas condiciones puede realizarse fácilmente mirando el gráfico de la serie temporal. Así, una serie temporal que exhiba una marcada tendencia creciente tendrá una media también creciente en el tiempo, por lo que lo más probable es que el proceso estocástico que ha generado dicha serie temporal no sea estacionario en media; del mismo modo, una serie temporal que muestre fluctuaciones de amplitud desigual en el tiempo seguramente no procederá de un proceso estocástico estacionario en varianza. La diferencia entre ambos tipos de series queda patente en los gráficos siguientes.

Pero, ¿por qué resulta importante para el investigador que el proceso analizado sea estacionario? La razón fundamental es que los modelos de predicción de series temporales que veremos a continuación están diseñados para ser utilizados con procesos de este tipo. Si las características del proceso cambian a lo largo del tiempo, resultará difícil representar la serie para intervalos de tiempo pasados y futuros mediante un modelo lineal sencillo, no pudiéndose por tanto realizar previsiones fiables para la variable en estudio.

Sin embargo, por regla general, las series económicas no son series que procedan de procesos estacionarios, sino que suelen tener una tendencia, ya sea creciente o decreciente, y variabilidad no constante. Dicha limitación en la práctica no es tan importante porque las series no estacionarias se pueden transformar en otras aproximadamente estacionarias después de aplicar diferencias a la serie en una ó más etapas. Por ello, cuando estemos analizando una serie económica que no sea estacionaria en media, deberemos trabajar con la serie en diferencias, especificando y estimando un modelo para la misma. Si además, observamos que la serie presenta no estacionariedad en varianza, deberemos transformarla tomando logaritmos antes de aplicar diferencias en la serie.

Las transformaciones de Box-Cox son una familia de transformaciones potenciales usadas en estadística para corregir sesgos en la distribución de errores, para corregir varianzas desiguales (para diferentes valores de la variable predictora) y principalmente para corregir la no linealidad en la relación (mejorar correlación entre las variables). La transformación de Box-Cox, gira en torno a al parámeto \(\lambda\):

\[\begin{equation}\tilde {y} _t = \left\lbrace\begin{array}{ll} log (y_t) & \lambda = 0 (y_t ^ \lambda - 1)\\ (y_t ^ \lambda - 1)\lambda & \lambda \neq 0 \\ \end{array}\right.\end{equation}\]

En R se realiza la transformación box-cox con la función “boxcox”" de la libreria “EnvStats”. Por defecto selecciona el \(\lambda\) en base al criterio de Probability Plot Correlation Coefficient (PPCC).

library(EnvStats)

## 
## Attaching package: 'EnvStats'

## The following object is masked from 'package:MASS':
## 
##     boxcox

## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm

## The following object is masked from 'package:base':
## 
##     print.default

y <- rlnormAlt(30, mean = 10, cv = 2)
boxcox(y, objective.name = "Log-Likelihood") 

## 
## Results of Box-Cox Transformation
## ---------------------------------
## 
## Objective Name:                  Log-Likelihood
## 
## Data:                            y
## 
## Sample Size:                     30
## 
##  lambda Log-Likelihood
##    -2.0     -164.01533
##    -1.5     -140.80593
##    -1.0     -120.59298
##    -0.5     -105.27508
##     0.0      -98.44032
##     0.5     -104.14089
##     1.0     -122.02228
##     1.5     -147.57348
##     2.0     -177.03427

Posteriormente, la predicción que realicemos con las series transformadas habrá que traducirla a una predicción para la serie original, en cuyo análisis estaba interesado inicialmente el investigador, deshaciendo las diferencias y aplicando antilogaritmos según convenga.

Por último, antes de continuar avanzando, debemos hacer mención a un tipo de proceso estacionario particular: es el denominado ruido blanco, un proceso estocástico en el que las variables aleatorias que lo forman no están correlacionadas entre sí, siendo su esperanza matemática igual a cero y su varianza constante e igual a \(\sigma^2\).

En particular, supondremos que los errores de los procesos que veremos a continuación son ruidos blancos gaussianos, formados por una sucesión de variables aleatorias con distribución Normal, esperanza cero, varianza constante e incorrelacionadas serialmente entre sí. Es decir, si \(e_t\) es ruido blanco gaussiano, entonces \(e_t \sim N(0,\sigma^2)\) y \(cov(e_t,e_{t+k})=0\).

En la figura siguiente se muestra la representación de un ruido blanco, en la que se puede apreciar claramente la estacionariedad de este proceso:

Operador de Retardos y Operador Diferencia

Antes de seguir avanzando, debemos mencionar dos operadores que utilizaremos frecuentemente a lo largo del capítulo. Por un lado, se define el operador de retardos, que denotaremos por B, como aquel operador que al ser aplicado a la serie la transforma de tal forma que: \(BY_t = Y_{t-1}\)

Es decir, el resultado de aplicar el operador \(B\) corresponde a retardar las observaciones un período.

Aplicada dos veces sobre la variable \(Y_t\) tendremos que: \(B(BY_t) = B^2Y_t = Y_{t-2}\)

y, en general, podemos decir que el operador \(B^k\) aplicado sobre una variable en el periodo t, la retarda k períodos tal que: \(B^kY_t = Y_{t-k}\)

Por su parte, el operador diferencia, el cual denotaremos por \(\Delta\), aplicado a una serie la transformará de tal forma que: \(\Delta Y_t = Y_t – Y_{t-1} = (1 – B) Y_t\)

Si aplicamos el operador diferencia dos veces a la serie tendremos que:

\[\Delta^2 Y_t =\Delta (\Delta Y_t) = \Delta (Y_t – Y_{t-1}) = \Delta Y_t – \Delta Y_{t-1} = Y_t – 2Y_{t-1} +Y{t-2} = (1–B)^2Y_t\]

Y en general, podemos escribir: \(\Delta^kY_t = (1 – B)^kY_t\)

Por lo que resulta evidente que la relación existente entre el operador diferencia y el operador retardo es: \(\Delta^k = (1 – B)^k\)

Modelización univariante de series temporales

La representación formal de los procesos aleatorios que generan series reales se puede realizar mediante modelos lineales de series temporales. Considerando que una determinada serie temporal ha sido generada por un proceso estocástico, en este epígrafe pasamos a describir los posibles modelos teóricos que permiten explicar el comportamiento de la misma y, por tanto, el de su proceso generador.

Las estructuras estocásticas estacionarias lineales que se tratarán de asociar a una serie de datos económicos se clasifican en tres tipos: modelos autorregresivos, modelos de medias móviles y modelos mixtos, los cuales pasamos a ver a continuación.

# Simulamos AR(1) con \phi_1=0.9
AR.1<-arima.sim(model=list(ar=c(.9)), n=100)
ts.plot(AR.1)

plot(acf(AR.1, type="correlation", plot=T))

plot(acf(AR.1, type="partial", plot=T))

# Simulamos AR(1) con \phi_1=-0.9
AR.1<-arima.sim(model=list(ar=c(-.9)), n=100)
ts.plot(AR.1)

plot(acf(AR.1, type="correlation", plot=T))

plot(acf(AR.1, type="partial",plot=T))

# Simulamos AR(2) con \phi_1=0.5 y \phi_2=0.2
AR.2<-arima.sim(model=list(ar=c(0.5,.2)),n=100)
ts.plot(AR.2)

plot(acf(AR.2,type="correlation",plot=T))

plot(acf(AR.2,type="partial",plot=T))

# Simulamos AR(2) con \phi_1=0.5 y \phi_2=-0.2
AR.2<-arima.sim(model=list(ar=c(.5,-.2)),n=100)
ts.plot(AR.2)

plot(acf(AR.2,type="correlation",plot=T))

plot(acf(AR.2,type="partial",plot=T))

# Simulamos MA(1) con \theta_1=-0.7
MA.1<-arima.sim(model=list(ma=c(-.7)),n=100)
ts.plot(MA.1)

plot(acf(MA.1,type="correlation",plot=T))

plot(acf(MA.1,type="partial",plot=T))

# Simulamos MA(1) con \theta_1=+0.7
MA.1<-arima.sim(model=list(ma=c(.7)),n=100)
ts.plot(MA.1)

plot(acf(MA.1,type="correlation",plot=T))

plot(acf(MA.1,type="partial",plot=T))

# Simulamos MA(2) con \theta_1=-0.7 y \theta_2=0.1
MA.2<-arima.sim(model=list(ma=c(-.7,.1)),n=100)
ts.plot(MA.2)

plot(acf(MA.2,type="correlation",plot=T))

plot(acf(MA.2,type="partial",plot=T))

# Simulamos MA(2) con \theta_1=+0.7 y \theta_2=0.1
MA.2<-arima.sim(model=list(ma=c(.7,.1)),n=100)
ts.plot(MA.2)

plot(acf(MA.2,type="correlation",plot=T))

plot(acf(MA.2,type="partial",plot=T))

# Simulamos MA(1) con \phi_1=0.7 y \theta_1=-0.1
ARMA.1<-arima.sim(model=list(ma=c(.7),ar=c(.1)),n=100)
ts.plot(ARMA.1)

plot(acf(ARMA.1,type="correlation",plot=T))

plot(acf(ARMA.1,type="partial",plot=T))

Fases para la elaboración de modelos univariantes

Aunque las fases que vamos a explicar a continuación son válidas para cualquier tipo de proceso univariante, vamos a centrarnos fundamentalmente en los procesos de tipo ARIMA. Básicamente, se trata de buscar un proceso ARIMA que de forma verosímil haya podido generar la serie temporal, es decir, que se adapte mejor a las características de la misma. Para ello seguiremos las siguientes fases:

1. Fase de identificación.

2. Fase de estimación

3. Fase de validación

4. Fase de predicción

Ejemplos prácticos

data(AirPassengers)
x<-AirPassengers

x =log(x)

par(mfcol = c(1, 2))
acf(x) 
pacf(x)

par(mfcol = c(1, 1))
x1 <- diff(x)
plot(x1, main="Pasajeros en Líneas Aéreas. Diferencias Log.")

par(mfcol = c(1, 2))
acf(x1) 
pacf(x1)

par(mfcol = c(1, 1))
x1_12 <- diff(x1, 12)
plot(x1_12, main="Pasajeros en Líneas Aéreas. Diferencias Log. Estacionales")

par(mfcol = c(1, 2))
acf(x1_12) 
pacf(x1_12)

library (forecast)
mod1 <- arima(x,c(1,1,0),c(0,1,1))
mod1

## 
## Call:
## arima(x = x, order = c(1, 1, 0), seasonal = c(0, 1, 1))
## 
## Coefficients:
##           ar1     sma1
##       -0.3395  -0.5619
## s.e.   0.0822   0.0748
## 
## sigma^2 estimated as 0.001367:  log likelihood = 243.74,  aic = -481.49

plot(mod1)

mod2 <- arima(x,c(0,1,1),c(0,1,1))
mod2

## 
## Call:
## arima(x = x, order = c(0, 1, 1), seasonal = c(0, 1, 1))
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 estimated as 0.001348:  log likelihood = 244.7,  aic = -483.4

plot (mod2)

mod3 <- arima(x,c(1,1,1),c(0,1,1))
mod3

## 
## Call:
## arima(x = x, order = c(1, 1, 1), seasonal = c(0, 1, 1))
## 
## Coefficients:
##          ar1      ma1     sma1
##       0.1960  -0.5784  -0.5643
## s.e.  0.2475   0.2132   0.0747
## 
## sigma^2 estimated as 0.001341:  log likelihood = 244.95,  aic = -481.9

plot(mod3)

par(mfcol = c(1, 2))
plot(x)
lines(x-mod1$residuals,col="red")
plot(mod1$residuals)

par(mfcol = c(1, 2))
plot(x)
lines(x-mod2$residuals,col="red")
plot(mod1$residuals)

par(mfcol = c(1, 2))
plot(x)
lines(x-mod3$residuals,col="red")
plot(mod1$residuals)

par(mfcol = c(1, 2))
acf(mod1$residuals)
pacf(mod1$residuals)

par(mfcol = c(1, 2))
acf(mod2$residuals)
pacf(mod2$residuals)

par(mfcol = c(1, 2))
acf(mod3$residuals)
pacf(mod3$residuals)

summary(mod1$residuals)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -0.1114640 -0.0232686 -0.0004461  0.0004500  0.0238061  0.1032393

par(mfcol = c(1, 1))
hist(mod1$residuals)

library(tseries)
jarque.bera.test(mod1$residuals)

## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Jarque Bera Test
## 
## Data:                            mod1$residuals
## 
## Test Statistic:                  X-squared = 3.607033
## 
## Test Statistic Parameter:        df = 2
## 
## P-value:                         0.1647186

summary(mod2$residuals)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -0.1185903 -0.0191628 -0.0000515  0.0005731  0.0223273  0.1085126

par(mfcol = c(1, 1))
hist(mod2$residuals)

library(tseries)
jarque.bera.test(mod2$residuals)

## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Jarque Bera Test
## 
## Data:                            mod2$residuals
## 
## Test Statistic:                  X-squared = 5.226523
## 
## Test Statistic Parameter:        df = 2
## 
## P-value:                         0.07329508

summary(mod3$residuals)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.245e-01 -2.005e-02  8.487e-05  6.215e-04  2.350e-02  1.138e-01

par(mfcol = c(1, 1))
hist(mod3$residuals)

library(tseries)
jarque.bera.test(mod3$residuals)

## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Jarque Bera Test
## 
## Data:                            mod3$residuals
## 
## Test Statistic:                  X-squared = 8.147297
## 
## Test Statistic Parameter:        df = 2
## 
## P-value:                         0.01701519

library(descomponer)
data(ipi)
x<-ts(ipi,c(2000,1),frequency=12)

x =log(x)

library(forecast)
mod1 <- auto.arima(x,test="pp")
mod1

## Series: x 
## ARIMA(1,0,2)(1,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1     ma2    sar1     sma1
##       0.9307  -0.7020  0.1629  0.3013  -0.7317
## s.e.  0.0408   0.0944  0.0862  0.1468   0.1194
## 
## sigma^2 estimated as 0.003462:  log likelihood=192.41
## AIC=-372.82   AICc=-372.17   BIC=-355.35

plot(mod1)

Introducción

Se dice que se ajusta el modelo paramétrico cuando se estiman sus parámetros a partir de un conjunto de observaciones que siguen dicho modelo, de manera que pueden hacerse predicciones de nuevos valores de \(Y\) conocido el valor de \(X\), y tener información precisa acerca de la incertidumbre asociada a la estimación y a la predicción.

Sin embargo, si el modelo paramétrico no es el adecuado al análisis de datos que estamos realizando, puede llevar a conclusiones que queden muy alejadas de la realidad, dado que el modelo paramétrico conlleva un grado de exactitud en las afirmaciones que de él se derivan y que son adecuadas siempre y cuando se cumplan los supuestos básicos sobre los que se apoya su construcción teórica. De hecho, los modelos paramétricos presentan una estructura teórica tan rígida que no pueden adaptarse a muchos conjuntos de datos de los que hoy día se disponen para el análisis económico.

La econometría no paramétrica aparece como consecuencia de intentos por solucionar problemas que existen en la econometría paramétrica como, por ejemplo, la consistencia entre los datos y los principios de maximización, homocedasticidad, o la necesidad de asumir una determinada relación, por lo general de forma lineal, entre las variables de interés. Esta preocupación llevó a una serie de investigadores a utilizar formas funcionales flexibles para aproximarse a relaciones desconocidas entre las variables. El plantear formas funcionales flexibles requiere el conocimiento del valor esperado de la variable \(Y\), condicional en las otras, \(X\). Esto conlleva la necesidad de estimar la función de densidad de \(Y\) condicional en \(X\).

La econometría no paramétrica no parte de supuestos sobre la distribución de probabilidad de las variables bajo estudio, sino que trata de estimar dicha distribución para encontrar la media condicional y los momentos de orden superior (por ejemplo, la varianza) de la variable de interés. Una de las desventajas de este método es, sin embargo, la necesidad de contar con muestras muy grandes si es que se desea estimar la función de relación entre ambas variables de manera precisa. Además el tamaño de la muestra debe aumentar considerablemente conforme aumenta el número de variables involucradas en la relación.

Los modelos de regresión paramétricos suponen que los datos observados provienen de variables aleatorias cuya distribución es conocida, salvo por la presencia de algunos parámetros cuyo valor se desconoce.

\[y_i= \beta_0 + \beta_1 x_i + e_i,\ \ i=1,2,…,n\]

con \(e_i\approx N(0;\sigma^2)\)

Este es un modelo estadístico con tres parámetros desconocidos: \(\beta_0\); \(\beta_1\) y \(\sigma^2\).

Una formulación general de un modelo de regresión paramétrico es la siguiente:

\[y_i=m(x_i;\theta)+e_i,\ \ i=1,2,…, n,\ \ \theta \in \Theta \in \Re\]

Donde \(m(x_i;\theta)\) es una función conocida de \(X\) y de \(\theta\), que es desconocido, \(e_i\) es una variable aleatoria idénticamente distribuida con \(E(e_i)=0\) y \(Var(e_i)=\sigma^2\). El modelo de regresión lineal simple sería un caso particular con \(\theta=(\beta_0, \beta_1)\) y \(m(x_i;\theta)=\beta_0 + \beta_1 x_i\).

Donde los valores de la variable explicativa son conocidos, por lo que se dice que el modelo tiene diseño fijo, y dado que la varianza de los errores es constante el modelo es Homocedástico.

Considerando una variable aleatoria bivariante con densidad conjunta \((Y,X)\), cabe definir la función de regresión \(f(y,x)\) como \(m(x,\theta)=E(Y/X=x)\), es decir, el valor esperado de \(Y\) cuando \(X\) toma el valor conocido \(x\). Entonces \(E(Y/X)=m(X,\theta))\), y definiendo \(e=Y-m(X,\theta)\), se tiene que: \(Y=m(X,\theta)+e\), \(E(e/X)=0\) y \(Var(e/X)=\sigma^2\).

Se supone que se observan ahora \(n\) pares de datos \((y_i,x_i),\ \ i=1…n\). Estos datos siguen el modelo de regresión no paramétrico: \(y_i=m(x_i;\theta)+e_i\).

Una vez establecido el modelo, el paso siguiente consiste en estimarlo (o ajustarlo) a partir de las n observaciones disponibles. Es decir, hay que construir un estimador de la función de regresión y un estimador de la varianza del error. Los procedimientos de estimación se conocen como métodos de suavizado.

El abanico de técnicas disponibles para estimar no paramétricamente la función de regresión es amplísimo e incluye, entre otras, las siguientes:

• Ajuste local de modelos paramétricos. Se basa en hacer varios (o incluso infinitos, desde un punto de vista teórico) ajustes paramétricos teniendo en cuenta únicamente los datos cercanos al punto donde se desea estimar la función.

• Suavizado mediante splines. Se plantea el problema de buscar la función \(\hat m(x)\) que minimiza la suma de los cuadrados de los errores (\(e_i=y_i-\hat m(x)\)) más un término que penaliza la falta de suavidad de las funciones (\(\hat m(x)\)) candidatas (en términos de la integral del cuadrado de su derivada segunda).

• Métodos basados en series ortogonales de funciones. Se elige una base ortonormal del espacio vectorial de funciones y se estiman los coeficientes del desarrollo en esa base de la función de regresión. Los ajustes por series de Fourier y mediante wavelets son los dos enfoques más utilizados.

• Técnicas de aprendizaje supervisado. Las redes neuronales, los k vecinos más cercanos y los árboles de regresión se usan habitualmente para estimar \(m(x)\).

Ajuste local (estimadores núcleo)

Los histogramas son siempre, por naturaleza, funciones discontinuas; sin embargo, en muchos casos es razonable suponer que la función de densidad de la variable que se está estimando es continua. En este sentido, los histogramas son estimadores insatisfactorios. Los histogramas tampoco son adecuados para estimar las modas, a lo sumo, pueden proporcionar “intervalos modales”, y al ser funciones constantes a trozos, su primera derivada es cero en casi todo punto, lo que les hace completamente inadecuados para estimar la derivada de la función de densidad.

Los estimadores de tipo núcleo (o kernel) fueron diseñados para superar estas dificultades. La idea original es bastante antigua y se remonta a los trabajos de Rosenblatt y Parzen en los años 50 y primeros 60. Los estimadores kernel son, sin duda, los más utilizados y mejor estudiados en la teoría no paramétrica.

Dada un m.a.s (\(X_1,X_2,...X_n\)) con densidad \(f\), estimamos dicha densidad en un punto \(t\) por medio del estimador:

\[\hat f(t)=\frac{1}{nh}\sum_{i=1}^n K(\frac{t-X_i}{h})\]

donde \(h\) es una sucesión de parámetros de suavizado, llamados ventanas, amplitudes o ancho de banda (windows, bandwidths) que deben tender a cero ”lentamente" (\(h \Rightarrow 0\), \(nh\Rightarrow \infty\)) para poder asegurar que \(\hat f\) tiende a la verdadera densidad \(f\) de las variables \(X_i\) y \(K\) es una función que cumple \(\int K=1\). Por ejemplo:

• Nucleo Gaussiano \(\frac{1}{\sqrt{2 \pi}}e^{- \frac {u^2} 2}\)

• Núcleo Epanechnikov \(\frac 3 4 (1-u^2)I_{\mid u \mid < 1}\)

donde \(I_{\mid u \mid < 1}\) es la función que vale 1 si \({\mid u \mid < 1}\) y cero si \({\mid u \mid \geq 1}\)

• Núcleo Triangular \((1-\mid u \mid)I_{\mid u \mid < 1}\)

• Núcleo Uniforme \(\frac 1 2 I_{\mid u \mid < 1}\)

• Núcleo Biweight \(\frac {15} {16} (1-u^2)I_{\mid u \mid < 1}\)

• Núcleo Triweight \(\frac {35} {32} (1-u^2)I_{\mid u \mid < 1}\)

Los programas informáticos eligen el ancho de ventana, \(h\), siguiendo criterios de optimización (error cuadrático medio).

En R la estimación de una función de densidad kernel se realiza con la función “density”. Con los datos del vector \(x\), hay que realizar el siguiente programa:

x <- c(2.1,2.6,1.9,4.5,0.7,4.6,5.4,2.9,5.4,0.2)
density(x, kernel="epanechnikov")

## 
## Call:
##  density.default(x = x, kernel = "epanechnikov")
## 
## Data: x (10 obs.);   Bandwidth 'bw' = 1.065
## 
##        x                  y          
##  Min.   :-2.99424   Min.   :0.00000  
##  1st Qu.:-0.09712   1st Qu.:0.02366  
##  Median : 2.80000   Median :0.09427  
##  Mean   : 2.80000   Mean   :0.08621  
##  3rd Qu.: 5.69712   3rd Qu.:0.15245  
##  Max.   : 8.59424   Max.   :0.16948

plot(density(x, kernel="epanechnikov"))

Los estimadores núcleo establecen que el peso de (\(X_i\),\(Y_i\)) en la estimación de \(m(x)\) es:

\[W_i(t,X_i)=\frac {\frac 1 h K(\frac{t-X_i}{h})}{\hat f(t)}\]

donde \(K(t)\) es una función de densidad simétrica (por ejemplo, la normal estándar) y \(\hat f(t)\) es un estimador kernel de la densidad como el definido en el apartado anterior.

\(W_i(t,X_i)\) es, para cada i, una función de ponderación que da “mayor importancia” a los valores \(X_i\) de la variable auxiliar que están cercanos a t.

Una expresión alternativa para \(W_i(t,X_i)\)

\[W_i(t,X_i)=\frac {K(\frac{t-X_i}{h})}{\sum_{j=1}^nK(\frac{t-X_i}{h})}\]

A partir de los pesos puede resolverse el problema de mínimos cuadrados ponderados siguiente:

\[min_{a,b}\sum_{i=1}^n W_i(Y_i-(\beta_0+\beta_1(t-X_i)))^2\]

los parámetros así obtenidos dependen de \(t\), porque los pesos \(W_i\) también dependen de \(t\). La recta de regresión localmente ajustada alrededor de \(t\) sería:

\[l_t(X)=\hat \beta_0(t)+\hat \beta_1(t)(t-X_i)\]

Y la estimación de la función en el punto en donde \(X=t\)

\[\hat m(t) =l_t(t)=\hat \beta_0(t)\]

Las funciones núcleo usadas en la estimación no paramétrica de la regresión son las mismas que en la densidad.

Si se generaliza al ajuste local de regresiones polinómicas de mayor grado, es decir si pretendemos estimar una forma lineal del tipo:

\[\beta_0+\beta_1 X_i+\beta_3 X_i^2+...+\beta_q X_i^q\]

con la salvedad de que en vez del valor \(X_i\) en la regresión lineal múltiple se utiliza el valor \(t-X_i\). El estimador de polinomios locales de grado \(q\), asignando los pesos obtenidos mediante la función núcleo. se resuelve mediante el siguiente problema de regresión polinómica ponderada:

\[min_{a,b}\sum_{i=1}^n W_i(Y_i-(\beta_0+\beta_1(t-X_i)+\beta_2(t-X_i)^2+...+\beta_q(t-X_i)^q))^2\]

Los parámetros \(\hat \beta_j= \hat \beta_j(t)\) dependen del punto t en donde se realiza la estimación, y el polinomio ajustado localmente alrededor de t sería:

\[P_{q,t}=\sum_{j=0}^q \hat \beta_j(t-X_i)^j\]

Siendo \(m(t)\) el valor de dicho polinomio estimado en el punto en donde \(X=t\):

\[\hat m(t) =P_{q,o}= \hat \beta_0(t)\]

En el caso particular del ajuste de un polinomio de grado cero, se obtiene el estimador de Nadaraya−Watson, o estimador núcleo de la regresión:

\[\hat m_k(t)=\frac {\sum_{i=0}^n K(\frac{t-X_i}{h})Y_i}{\sum_{i=0}^n K(\frac{t-X_i}{h})}\]

El estimador del parámetro de suavizado \(h\) tiene una importancia crucial en el aspecto y propiedades del estimador de la función de regresión. Valores pequeños de \(h\) dan mayor flexibilidad al estimador y le permiten acercarse a todos los datos observados, pero originan altos errores de predicción (sobre-estimación). Valores mas altos de \(h\) ofrecerán un menor grado de ajuste a los datos pero predicen mejor, aunque si \(h\) es demasiado elevado tendremos una falta de ajuste a los datos (sub-estimación).

En la figura siguiente se realizar la representación gráfica de la regresión kernel realizada con el estimador de Nadaraya–Watson con diferentes parámetros de suavizado a la serie temporal co2.

plot (time(co2),co2, main="Regresión kernel con diferentes parámetros de suavizado",type="l")
lines(ksmooth(time(co2),co2, "normal", bandwidth = 2), col = 2)
lines(ksmooth(time(co2),co2, "normal", bandwidth = 5), col = 3)

La función “dpill” de la libreria “KernSmooth”, permite una selección automática del bandwidth más idoneo, siempre que no existan fuertes irregularidades en la serie de datos originales (outliers).

library(KernSmooth)

## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009

plot(time(co2),co2, main="Regresión kernel con función dpill",type="l")
lines(ksmooth(time(co2),co2, "normal", bandwidth = dpill(cars$speed, cars$dist)), col = 2)

Regresión LOESS

La Regresión LOESS (inicialmente Lowess: Locally Weighted Scatterplot Smoothing) es un modelo de regresión no-paramétrica que también se conoce como regresión local.

LOESS combina la sencillez de la regresión lineal por mínimos cuadrados con la flexibilidad de la regresión no lineal. Mediante el ajuste de modelos sencillos sobre subconjuntos locales de datos, crea una función que describe la parte determinista de la variación en los datos punto a punto. Uno de los principales atractivos de este método es que no resulta necesario especificar una función global para ajustar un modelo a los datos.

La regression LOESS, fue propuesta originalmente por Cleveland (1979) y desarrollada por Cleveland y Devlin (1988).

Dado un conjunto de pares de datos, \((x_t,y_t)\), correspondientes a observaciones de las variables \((X,Y)\). El objetivo de un suavizado en un diagrama de dispersión es estimar la componente determinística,\(\mu(\cdot)\) cuando se utiliza el siguiente modelo de la relación entre las variables:

\[y_t=\mu(x_t)+e_t\]

Un modelo de regresió lineal simple, supone estimar:

\[\mu(X_t)=\beta_0+\beta_1x_t\]

con todas las observaciones participando en la estimación.

La regresión local es un enfoque de ajuste de curvas (o superficies) a datos mediante suavizados en los que el ajuste en \(x\) se realiza utilizando únicamente observaciones en un entorno de \(x\). Al realizar una regresión local puede utilizarse una familia paramétrica al igual que en un ajuste de regresión global pero solamente se realiza el ajuste localmente.

En la práctica se realizan ciertas suposiciones:

• sobre la función de regresión \(\mu(\cdot)\) tales como continuidad y derivabilidad de manera que pueda estar bien aproximada localmente por polinomios de un cierto grado.

• sobre la variabilidad de \(Y\) alrededor de la curva \(\mu(\cdot)\).

Los métodos de estimación que resultan de este tipo de modelos, requieren que para cada punto \(x\), se defina un entorno. Dentro de ese entorno suponemos que la función regresora es aproximada por una función paramétrica polinómica.Por ejemplo, un polinomio cuadrático: \(g(u) = a_0 + a_1 (u - x) +a_2 (u -x)^2\).

Se estiman los parámetros con las observaciones en ese entorno, denominándose ajuste local a la función ajustada evaluada en \(x\).

Generalmente, se incorpora una función de peso, \(w(u)\), para dar mayor peso a los valores \(x_t\) que se encuentran cerca de \(x\). Los criterios de estimación dependen de los supuestos que se realicen sobre la distribución de \(Y\) dado \(x\). Si, por ejemplo, suponemos que tiene distribución Gaussiana con varianza constante, es adecuado utilizar el Método de Mínimos Cuadrados.

Los subconjuntos de los datos utilizados para el ajuste por mínimos cuadrados ponderados acaban estando determinados por un parámetro de suavización que define el ancho de banda, \(\alpha\). Este parámetro es un número entre \(\frac {(\lambda + 1)}{n}\) y 1, donde \(\lambda\) denota el grado del polinomio local. El valor de \(\alpha\) es la proporción de los datos utilizados en cada ajuste. A \(\alpha\) se le llama parámetro de suavización porque controla la flexibilidad de la función de regresión. Valores grandes de \(\alpha\) producen curvas suaves; valores pequeños hacen que la curva se ajuste tal vez demasiado a los datos. En ocasiones se recomienda utilizar valores en el rango que va de \(0.25\) a \(0.5\).

Este es el codigo en R para el suavizado con Loess.

co2.l=loess(co2~time(co2))
summary(co2.l)

## Call:
## loess(formula = co2 ~ time(co2))
## 
## Number of Observations: 468 
## Equivalent Number of Parameters: 4.34 
## Residual Standard Error: 2.113 
## Trace of smoother matrix: 4.72  (exact)
## 
## Control settings:
##   span     :  0.75 
##   degree   :  2 
##   family   :  gaussian
##   surface  :  interpolate      cell = 0.2
##   normalize:  TRUE
##  parametric:  FALSE
## drop.square:  FALSE

plot (co2, main="Regresión LOESS")
lines (ts(co2.l$fitted, frequency = 12, start = c(1959, 1)), col=2)

Splines

Introducción

Para poder estimar la función \(f\) de la forma más sencilla posible, deberíamos poder representarla de forma que \(Y_i=f(x_i)+e_i\), \(i=1,2,...,n\) se convierta en un modelo lineal.

Y esto se puede hacer eligiendo una base de funciones de dimensión \(q\) que genere un subespacio de funciones que incluya a \(f\) como elemento y que pueda expresarse como:

\[f(x)=\sum_{j=1}^q \beta_j s_j(x)\]

siendo \(\beta_j\) un parámetro desconocido, asociado al elemento \(j\), \(s_j(x)\) de dicha base de funciones, de manera que:

\[Y_i=\sum_{j=1}^q \beta_j s_j(x) + e_i\] (1)

se convierte en un modelo lineal de dimensión \(q\).

La regresión con funciones base polinómicas es la propuesta más sencilla para este tipo de estimaciones.

Supongamos que \(f\) es un polinomio de grado 4 de forma que el espacio de polinomios de grado 4 contiene a \(f\). Una base de este subespacio es:

\[\begin{equation} \left\lbrace\begin{array}{ll} s_1(x)=1\\ s_2(x)=x\\ s_3(x)=x^2\\ s_4(x)=x^3\\ s_5(x)=x^4 \end{array}\right.\end{equation}\]

De manera que el modelo (1) se convierte en:

\[Y_i=\beta_1 + \beta_2 x_i + \beta_3 x_i^2 + \beta_4 x_i^3 + \beta_5 x_i^4 + e_i\]

Un spline es una curva diferenciable definida en porciones mediante polinomios, que se utiliza como base de funciones para aproximar curvas con formas complicadas.

Una función spline está formada por varios polinomios, cada uno definido sobre un subintervalo, que se unen entre sí obedeciendo a ciertas condiciones de continuidad.

Supongamos que se ha fijado un entero \(q \neq 0\), de manera que disponemos de \(q+1\) puntos, a los que denominaremos nodos (knots), tales que \(t_0<t_1<t_2<...<t_q\), en los que troceamos nuestro conjunto de datos. Decimos entonces que una función spline de grado \(q\) con nodos en \(t_0,t_1,t_2,...,t_q\) es una función \(S\) que satisface las condiciones:

1. en cada intervalo \([t_{j-1}, t_j)\), \(S\) es un polinomio de grado menor o igual a \(q\).

2. \(S\) tiene una derivada de orden \((q-1)\) continua en \([t_0,t_q]\).

La flexibilidad de estos modelos se puede controlar de dos formas:

1. Con el número de puntos de corte (knots) introducidos. A mayor número mayor flexibilidad.

2. Con el grado del polinomio empleado.

Los splines de grado \(0\) son funciones constantes por zonas. La expresión matemática de un spline de grado \(0\) es la siguiente:

\[\begin{equation} \left\lbrace\begin{array}{ll} s_0(x)=c_0 & x \in [t_0,t_1) \\ s_j(x)=c_j & x \in [t_{j-1},t_j) \\ ...\\ s_{q-1}(x)=c_{q-1} & x \in [t_{q-1},t_q) \end{array}\right.\end{equation}\]

En la figura siguiente se muestran las gráficas correspondientes a los splines de grado cero:

Los splines de grado \(0\), se definen en un solo tramo de nodo, y ni siquiera es continua en los nodos. Equivale a realizar una regresión por tramos.

\[Y_i=\beta_0 c_0 x_i + \beta_1 c_1 x_i + ... + \beta_{q-1} c_{q-1} x_i+ e_i\]

donde \(c_j\) toma valor \(1\) en \([t_{j-1},t_j)\) y cero en el resto.

Un spline de grado 1 o lineal se expresa como:

\[\begin{equation} \left\lbrace\begin{array}{ll} s_0(x)=a_0x + b_0 & x \in [t_0,t_1) \\ s_j(x)=a_jx + b_j & x \in [t_{j-1},t_j) \\ ...\\ s_{q-1}(x)=a_{q-1}x + b_{q-1} & x \in [t_{q-1},t_q) \end{array}\right.\end{equation}\]

La representación gráfica de un spline lineal aparece en la figura siguiente:

Las funciones de spilines más comúnmente utilizadas son las de grado tres ó cúbicas. Son polinomios de grado tres a trozos, que son continuos en los nodos al igual que su primera y segunda derivada, proporcionando un excelente ajuste a los puntos tabulados y a través de un cálculo que no es excesivamente complejo.

Sobre cada intervalo \([t_{j-1},t_j)\), \(S\) está definido por un polinomio cúbico diferente.

\[\begin{equation} \left\lbrace\begin{array}{ll} s_0(x)=a_0x^3 + b_0x^2 + c_0x + d_0 & x \in [t_0,t_1) \\ s_j(x)=a_jx^3 + b_jx^2 + c_jx + d_j & x \in [t_{j-1},t_j) \\ ...\\ s_{q-1}(x)=a_{q-1}x^3 + b_{q-1}x^2 + c_{q-1}x + d_{q-1} & x \in [t_{q-1},t_q) \end{array}\right.\end{equation}\]

Los polinomios \(S_{j-1}\) y \(S_j\) interpolan el mismo valor en el punto \(t_j\), para que se cumpla: \(S_{j-1}(x_i)=y_i=S_j(x_i)\), por lo que se garantiza que \(S\) es continuo en todo el intervalo. Además, se supone que \(S'\) y \(S''\) son continuas, condición que se emplea en la deducción de una expresión para la función del spline cúbico.

Aplicando las condiciones de continuidad del spline \(S\) y de las derivadas primera \(S'\) y segunda \(S''\), es posible encontrar la expresión analítica del spline.

Los elementos de una base de splines cúbicos son polinomios de grado 3. Un Spline cúbico se representa en la figura siguiente:

# LEEMOS DATOS
library(ISLR)
data("Wage")
head(Wage)

##        year age           maritl     race       education
## 231655 2006  18 1. Never Married 1. White    1. < HS Grad
## 86582  2004  24 1. Never Married 1. White 4. College Grad
## 161300 2003  45       2. Married 1. White 3. Some College
## 155159 2003  43       2. Married 3. Asian 4. College Grad
## 11443  2005  50      4. Divorced 1. White      2. HS Grad
## 376662 2008  54       2. Married 1. White 4. College Grad
##                    region       jobclass         health health_ins
## 231655 2. Middle Atlantic  1. Industrial      1. <=Good      2. No
## 86582  2. Middle Atlantic 2. Information 2. >=Very Good      2. No
## 161300 2. Middle Atlantic  1. Industrial      1. <=Good     1. Yes
## 155159 2. Middle Atlantic 2. Information 2. >=Very Good     1. Yes
## 11443  2. Middle Atlantic 2. Information      1. <=Good     1. Yes
## 376662 2. Middle Atlantic 2. Information 2. >=Very Good     1. Yes
##         logwage      wage
## 231655 4.318063  75.04315
## 86582  4.255273  70.47602
## 161300 4.875061 130.98218
## 155159 5.041393 154.68529
## 11443  4.318063  75.04315
## 376662 4.845098 127.11574

# AUSTE DEL MODELO
library(splines)
modelo_splines <- lm(wage ~ bs(age, knots = c(25,40,60), degree = 3),data = Wage)
summary(modelo_splines)

## 
## Call:
## lm(formula = wage ~ bs(age, knots = c(25, 40, 60), degree = 3), 
##     data = Wage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -98.832 -24.537  -5.049  15.209 203.207 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                   60.494      9.460   6.394
## bs(age, knots = c(25, 40, 60), degree = 3)1    3.980     12.538   0.317
## bs(age, knots = c(25, 40, 60), degree = 3)2   44.631      9.626   4.636
## bs(age, knots = c(25, 40, 60), degree = 3)3   62.839     10.755   5.843
## bs(age, knots = c(25, 40, 60), degree = 3)4   55.991     10.706   5.230
## bs(age, knots = c(25, 40, 60), degree = 3)5   50.688     14.402   3.520
## bs(age, knots = c(25, 40, 60), degree = 3)6   16.606     19.126   0.868
##                                             Pr(>|t|)    
## (Intercept)                                 1.86e-10 ***
## bs(age, knots = c(25, 40, 60), degree = 3)1 0.750899    
## bs(age, knots = c(25, 40, 60), degree = 3)2 3.70e-06 ***
## bs(age, knots = c(25, 40, 60), degree = 3)3 5.69e-09 ***
## bs(age, knots = c(25, 40, 60), degree = 3)4 1.81e-07 ***
## bs(age, knots = c(25, 40, 60), degree = 3)5 0.000439 ***
## bs(age, knots = c(25, 40, 60), degree = 3)6 0.385338    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.92 on 2993 degrees of freedom
## Multiple R-squared:  0.08642,    Adjusted R-squared:  0.08459 
## F-statistic: 47.19 on 6 and 2993 DF,  p-value: < 2.2e-16

# INTERPOLACIÓN DE PUNTOS DENTRO DEL RANGO DEL PREDICTOR
limites <- range(Wage$age)
nuevos_puntos <- seq(from = limites[1], to = limites[2], by = 1)
nuevos_puntos <- data.frame(age = nuevos_puntos)
head(nuevos_puntos)

##   age
## 1  18
## 2  19
## 3  20
## 4  21
## 5  22
## 6  23

tail(nuevos_puntos)

##    age
## 58  75
## 59  76
## 60  77
## 61  78
## 62  79
## 63  80

# PREDICCIÓN DE LA VARIABLE RESPUESTA Y DEL ERROR ESTÁNDAR
predicciones <- predict(modelo_splines, newdata = nuevos_puntos, se.fit = TRUE,level = 0.95)

# CÁLCULO DEL INTERVALO DE CONFIANZA SUPERIOR E INFERIOR
intervalo_conf <- data.frame(inferior = predicciones$fit - 1.96*predicciones$se.fit,superior = predicciones$fit +1.96*predicciones$se.fit)
attach(Wage)
plot(x = age, y = wage, pch = 20, col = "darkgrey")
title("Cubic Spline, knots: 25, 40, 60")
lines(x = nuevos_puntos$age, predicciones$fit, col = "red", lwd = 2)
lines(x = nuevos_puntos$age, intervalo_conf$inferior, col = "blue",lwd = 2, lty = 2)
lines(x = nuevos_puntos$age, intervalo_conf$superior, col = "blue",lwd = 2, lty = 2)

# LEEMOS DATOS
library(descomponer)
data(ipi)
x<-ts(ipi,c(2000,1),frequency=12)

# CUBIC SPLINE CON 3 KNOTS EQUIDISTRIBUIDOS 
ipi_splines <- lm(x ~ bs(time(x), df = 6, degree = 3))
summary(ipi_splines)

## 
## Call:
## lm(formula = x ~ bs(time(x), df = 6, degree = 3))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.812  -5.573   1.883   7.147  20.032 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        99.595      5.682  17.529   <2e-16 ***
## bs(time(x), df = 6, degree = 3)1  -10.610     10.653  -0.996   0.3210    
## bs(time(x), df = 6, degree = 3)2   13.940      6.992   1.994   0.0481 *  
## bs(time(x), df = 6, degree = 3)3   10.999      8.511   1.292   0.1983    
## bs(time(x), df = 6, degree = 3)4   -9.476      7.926  -1.195   0.2339    
## bs(time(x), df = 6, degree = 3)5    5.248      8.887   0.591   0.5558    
## bs(time(x), df = 6, degree = 3)6  -10.963      7.913  -1.385   0.1681    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.39 on 141 degrees of freedom
## Multiple R-squared:  0.2009, Adjusted R-squared:  0.1669 
## F-statistic: 5.907 on 6 and 141 DF,  p-value: 1.594e-05

# GRÁFICA DEL SUAVIZADO
x.spline<-ts(ipi_splines$fitted.values,c(2000,1),frequency=12)
plot(x,main = "IPI de Cantabria")
lines(x.spline,col="red")

attach(cars)
plot(speed, dist, main = "data(cars) & smoothing splines")
cars.spl1 <- smooth.spline(speed, dist)
cars.spl1

## Call:
## smooth.spline(x = speed, y = dist)
## 
## Smoothing Parameter  spar= 0.7801305  lambda= 0.1112206 (11 iterations)
## Equivalent Degrees of Freedom (Df): 2.635278
## Penalized Criterion (RSS): 4187.776
## GCV: 244.1044

cars.spl2 <- smooth.spline(speed, dist, spar=0.10)
lines(cars.spl1, col = "blue")
lines(cars.spl2, col = "red")

plot(x, main = "IPI de Cantabria")
IPI.spl1 <- smooth.spline(x)
IPI.spl1

## Call:
## smooth.spline(x = x)
## 
## Smoothing Parameter  spar= 0.7554373  lambda= 0.0007791226 (11 iterations)
## Equivalent Degrees of Freedom (Df): 8.414817
## Penalized Criterion (RSS): 16314.77
## GCV: 123.9264

IPI.spl2 <- smooth.spline(x, spar=0.50)
lines(IPI.spl1, col = "blue")
lines(IPI.spl2, col = "red")

Aproximación por series de Fourier

La forma de Fourier permite aproximar arbitrariamente cerca, tanto a la función como a sus derivadas, sobre todo el dominio de definición de las mismas. La idea que subyace en este tipo de aproximaciones (que podrían denominarse semi-no-paramétricas) es ampliar el orden de la base de expansión, cuando el tamaño de la muestra aumenta, hasta conseguir la convergencia asintótica de la función aproximante a la verdadera función generadora de los datos y a sus derivadas (Gallant, A.R.; 1981, 1984).

Un polinomio de Fourier viene dado por la expresión:

\[\frac \eta 2 + \sum_{j=1}^k[a_j\cos(jw_0t)+b_j\sin(jw_0t)]\]

donde \(k\) es es el número de ciclos teóricos o armónicos que consideramos, siendo el máximo \(N/2\), \(w_0=\frac {2\pi} N\) es la frecuencia fundamental (también denominada frecuencia angular fundamental) y \(t\) toma los valores enteros comprendidos entre \(1\) y \(N\) (es decir, \(t = 1, 2, 3, ...N\)).

Los coeficientes de los armónicos vienen dados por las expresiones:

\[\frac \eta 2 = \frac 2 N \sum_{i=1}^N y_i\]

\[a_j=\frac 2 N \sum_{i=1}^N y_i \cos(jw_0t)\]

\[b_j=\frac 2 N \sum_{i=1}^N y_i \sin(jw_0t)\]

La aproximación a una función no periódica \(m(x_i;\theta)\) por una serie de expansión de Fourier se escribe como:

\[m(x_i;\theta)= \eta + \sum_{j=1}^J[a_j\cos(jx_i)+b_j\sin(jx_i)]\]

El vector de parámetros es \(\theta=(\eta, a_1, b_1,...,a_J,b_J)\) de longitud \(K=1+2J\).

Suponiendo que los datos siguieran el modelo \(y_i=m(x_i;\theta)+e_i\), para \(i=1,2,…,N\), estimariamos \(\theta\) por mínimos cuadrados minimizando:

\[S_n(\theta)=\frac 1 N \sum_{i=1}^N[y_i-m_K(x_i;\theta)]^2\]

Dado que la variable exógena \(x_i\) no esta expresada en forma periódica, debe de transformase o normalizarse en un intervalo de longitud menor que \(2\pi\), \([0,2\pi]\).

Como ejemplo, vamos a utilizar la base de datos de la Agencia Española de Meteorológica (Aemet) desde el R-package fda.usc. La base de datos contiene mediciones diarias de temperatura, velocidad del viento y precipitaciones de 73 diferentes estaciones meteorológicas de España para los años 1980 a 2009. En este ejemplo, vamos a analizar las temperaturas medias diarias de Santander, que representamos gráficamente en R con la siguiente programación:

library(fda) 

## Loading required package: Matrix

## 
## Attaching package: 'fda'

## The following object is masked from 'package:forecast':
## 
##     fourier

## The following object is masked from 'package:graphics':
## 
##     matplot

library(fda.usc)

## Loading required package: mgcv

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:forecast':
## 
##     getResponse

## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

## Loading required package: rpart

data(aemet, package = "fda.usc")
tt = aemet$temp$argvals
temp = as.data.frame(aemet$temp$data, row.names=F)
range.tt = aemet$temp$rangeval
inv.temp = data.frame(t(aemet$temp$data)) # 365 x 73 matrix
names(inv.temp) = aemet$df$name 
plot(ts(inv.temp[,21]),main="Temperaturas medias diarias Santander 1980-2009")

Para conocer los ciclos que dominan en una serie temporal hay que representar el peridograma. El periodograma mide la aportación de los componentes periodicos de las frecuencias de la serie temporal (armonicos) a su varianza total. Si el periodograma presenta un pico en una frecuencia determinada, indica que dicha frecuencia tiene más importancia en la dinámica de la serie que el resto. Utilizamos la función gperiodograma de la librería descomponer para representar el periodograma de

# GENERAMOS UNA SERIE TEMPORAL ARMONICA CON TENDENCIA, DOS CICLOS DE FRECUENCIA 2 y 4 Y ERROR ALEATORIO

t=seq(1:100)
x=0.5*t+5*cos(2*pi*t/50)+0.5*sin(2*pi*t/50)+5*cos(2*pi*t/25)+0.5*sin(2*pi*t/50)+rnorm(100,0,2)
plot(x,type="l")

# PERIODOGRAMA DE LA SERIE TEORICA

library(descomponer)
gperiodograma(x)

# PERIODOGRAMA DE LA SERIE DE TEMPERATURAS DE SANTANDER
gperiodograma(inv.temp[,21])

El periodograma de la serie de temperaturas medias diarias muestra que la serie está dominada por los ciclos de baja frecuencia.

A continuación se van a suavizar estas temperaturas diarias utilizando funciones periódicas de Fourier, en concreto vamos a utilizar las funciones de base igual a 5. Es decir, los armónicos que se obtendrían con:

\[\sum_{j=1}^5[a_j\cos(jw_0t)+b_j\sin(jw_0t)]\]

Santander5 = create.fourier.basis(rangeval = range(tt), nbasis = 5)
plot(Santander5)

Con la función smooth.basis(argvals=1:n, y, fdParobj) del R-package fda obtenemos la serie suavizada. Respecto a los parámetros, argvals es el dominio, y es el conjunto de valores a suavizar y fdParobj la función base utilizada como regresores:

Santanderfourier5.fd = smooth.basis(argvals = tt, y = inv.temp[,21], fdParobj = Santander5)
plot(ts(inv.temp[,21]), main="Temperaturas medias diarias Santander 1980-2009")
lines(Santanderfourier5.fd,col="red")

En el objeto que genera la función, incluye un data set “y2cMap” en donde se guardan las series de senos y cosenos que se utilizan en el suavizado, en “fd$coefs” estan los coeficientes, para obtener el resultdo del suavizado hay que multiplicar los coeficientes por los regresores. En el siguiente proceso se realiza la misma representación obteniendo el resultado del suavizado, que se ha transformado a objeto ts.

# Regresores
head(t(Santanderfourier5.fd$y2cMap))

##           const         sin1       cos1        sin2       cos2
## [1,] 0.05170402 0.0006310771 0.07311781 0.001262107 0.07310964
## [2,] 0.05170423 0.0019103806 0.07309602 0.003819475 0.07302161
## [3,] 0.05170487 0.0031891174 0.07305275 0.006372281 0.07284597
## [4,] 0.05170592 0.0044669066 0.07298802 0.008917483 0.07258292
## [5,] 0.05170740 0.0057433674 0.07290184 0.011452047 0.07223277
## [6,] 0.05170930 0.0070181196 0.07279423 0.013972954 0.07179594

# Coeficientes
Santanderfourier5.fd$fd$coefs

##             [,1]
## const 276.441022
## sin1  -38.390736
## cos1  -61.224525
## sin2    8.514355
## cos2   -1.642055

# Obtenermos resultado
resultado=(t(Santanderfourier5.fd$fd$coefs)%*%Santanderfourier5.fd$y2cMap)
head(ts(resultado[1,]))

## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1] 9.682978 9.657176 9.632933 9.610235 9.589071 9.569427

plot(ts(inv.temp[,21]),main="Temperaturas medias diarias Santander 1980-2009")
lines(ts(resultado[1,]),col="red")

Regresión band spectrum

Nerlove (1964) y Granger (1969) fueron los primeros investigadores en aplicar el analisis espectral a las series de tiempo en economía. El uso del analisis espectral requiere un cambio en el modo de ver las series económicas, al pasar de la perspectiva del tiempo al dominio de la frecuencia. El analisis espectral parte de la supsición de que cuanquier serie \(X_t\), puede ser transformada en ciclos formados con senos u cosenos:

\[X_t=\eta+\sum_{j=1}^N\left[a_j\cos \left(2\pi\frac{ft}n\right) +b_j\sin\left(2\pi\frac{ft}n\right)\right] \ \ \ \ \ \ (1)\]

donde \(\eta\) es la media de la serie, \(a_j\) y \(b_j\) son su amplitud, \(f\) son las frecuencias del conjunto de las \(n\) observaciones, \(t\) es un indice de tiempo que va de 1 a N, siendo N el numero de periodos para los cuales tenemos observaciones en el conjunto de datos.

El cociente \((\frac{ft}n)\) convierte cada valor de \(t\) en escala de tiempo en proporciones de \(2n\) y rango \(j\), desde \(1\) hasta \(n\), siendo \(n=\frac{N}2\) (es decir, 0,5 ciclos por intervalo de tiempo).

La dinámica de las altas frecuencias (los valores más altos de f) corresponden a los ciclos cortos en tanto que la dinámica de la bajas frecuencias (pequeños valores de f) van a corresponder con los ciclos largos. Si nosotros hacemos que \(\frac{ft}n=w\), la ecuación (1) quedaría así:

\[X_t=\eta+\sum_{j=1}^N[a_j\cos(\omega_j)+b_j\sin(\omega_j)] \ \ \ \ \ \ (2)\]

El análisis espectral puede utilizarse para identificar y cuantificar, en procesos aparentemente aperiodicos, sucesiones de ciclos de corto y largo plazo. Una serie dada \({X_t}\) puede contener diversos ciclos de diferentes frecuencias y amplitudes, y esa combinación de frecuencias y amplitudes de carácter cíclico la hacen aparecer como una serie no periódica e irregular. De hecho, la ecuación (2) muestra que cada observación \(t\) de una serie de tiempo, es el resultado de sumar los valores en \(t\) que resultan de \(N\) ciclos de diferente longitud y amplitud, a los que habría que añadir, si cabe, un término de error.

Una manera practica de pasar desde el dominio del tiempo al dominio de la frecuencia es premultiplicando los datos originales por una matriz ortogonal, \(W\), sugerida por Harvey (1978), con el elemento (j,t)th:

\[\begin{equation} w_{jt} = \left\lbrace\begin{array}{ll}\left(\frac{1}T\right) ^\frac{1}2 & \forall j=1\\ \left(\frac{2}T\right) ^\frac{1}2 \cos\left[\frac{\pi j(t-1)}T\right] & \forall j=2,4,6,..\frac{(T-2)}{(T-1)}\\ \left(\frac{2}T\right) ^\frac{1}2 \sin\left[\frac{\pi (j-1)(t-1)}T\right] & \forall j=3,5,7,..\frac{(T-2)}T\\ \left(\frac{1}T\right) ^\frac{1}2 (-1)^{t+1} & \forall j=T\end{array}\right.\end{equation}\]

La matriz \(W\) tiene la ventaja de ser ortogonal por lo que \(WW^T=I\).

Sea \(a\) un vector n x 1, el modelo transformado en el dominio de la frecuencia esta dado por: \(\hat a= Wa\).

Engle (1974), demostró que una regresión realizada con las series transformadas en el dominio de la frecuencia, Regresión Band Spectrum (RBS), no alteraba los supuestos básicos de la regresión clásica, cuyos estimadores eran Estimadores Lineales, Insesgados y Optimos (ELIO).

\[y=X\beta+u \ \ \ \ \ \ (3)\]

donde \(X\) es una matriz \(n x k\) de observaciones de \(k\) variables independientes, \(\beta\) es un vector \(k x I\) de parámetros, \(y\) es un vector \(n x 1\) de observaciones de la variable dependiente, y \(u\) en un vector \(n x I\) de pertubacciones de media cero y varianza constante, \(\sigma^2\).

El modelo se expresaría en el dominio de la frecuencia aplicando una transformación lineal a las variables dependiente e independientes, por ejemplo, premultiplicando todas las variables por la matriz ortogonal \(W\). La técnica de la RBS consiste en realizar el analisis de regresión en el dominio de la frecuencia, pero omitinedo determinadas oscilaciones periodicas. Con este procedimiento pueden tratarse problemas derivados de la estacionalidad de las series o de la autocorrelación en los residuos. Engle (1974) muestra que si los residuos están correlacionados serialmente y son generados por un procieso estacionario estocástico, la regresión en el dominio de la frecuencia es el estimador asintóticamente más eficiente de \(\beta\).

La transformación de la ecuación (3) del dominio del tiempo al dominio de la frecuencia en Engle (1974), se basa en la matriz \(W\), cuyo elemento \((t, s)\) esta dado por:

\[w_{ts}=\frac{1}{\sqrt n} e^{i\lambda_t s} \ \ \ \ \ \ s= 0,1,...,n-1\]

donde \(\lambda_t = 2\pi \frac{t}n\), \(t=0,1,...,n-1\), y \(i=\sqrt{-1}\).

Premultiplicando las observaciones de (3) por \(W\), obtenemos:

\[\dot y=\dot X\beta+\dot u \ \ \ \ \ \ (4)\]

donde \(\dot y = Wy\), \(\dot X = WX\),ç y \(\dot u = Wu\).

Si el vector de las perturbaciones en (3) cumple las hipótesis clásicas del modelo de regresión: \(E[u] = 0\) y \(E[uu']=\sigma^2 I_n\), entonces el vector de perturbaciones transformado al dominio de la frecuencia, \(\dot u\), tendrá las mimas propiedades.

Por otro lado, dado que la matriz W es ortogonal, \(WW^{T}= I\), entonces \(W^T\) sería la transpuesta de la completa conjugada de \(W\), de forma que las observaciones del modelo (4) acaban conteniendo el mismo tipo de información que las observaciones del modelo inicialmente planteado.

\[var(\dot u)=E(\dot u \dot u^T)=E(Wuu'W^T)=WE(uu')W^T=\sigma^2W \Omega W^T\]

si \(\Omega=I\) entonces \(var(\dot u)=\sigma^2I\).

Asuminendo que \(x\) es independiente de \(u\), el toerema de Gauss-Markov implicaría que:

\[\hat \beta = (\dot x' \dot x)^{-1}\dot x'\dot y\]

es el mejor estimador lineal insesgando (ELIO) de \(\dot \beta\). El estimador obtenido sería de hecho idéntico al estimador MCO de (3).

Estimar (4) manteniendo únicamente determinadas frecuencias, puede llevarse a cabo omitiendo las observaciones correspondientes a las restantes frecuencias, si bien, dado que las variables en (4) son complejas, Engle (1974) sugiere la transformada inversa de Fourier para recomponer el modelo estimado en términos de tiempo.

Utilizando el IPI de Cantabria se comprueba que los estimadores MCO coinciden con los estimadores de la regresión en el dominio de la frecuencia:

# DATOS IPI CANTABRIA
library(descomponer)
data("ipi")
t=seq(1:length(ipi))
t2=t^2

# REGRESIÓN MINIMO CUADRADA
summary(lm(ipi~t+t2))

## 
## Call:
## lm(formula = ipi ~ t + t2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.468  -4.931   2.215   6.942  23.111 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 95.2675702  2.9250060  32.570  < 2e-16 ***
## t            0.3413570  0.0906335   3.766  0.00024 ***
## t2          -0.0025658  0.0005892  -4.355  2.5e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.7 on 145 degrees of freedom
## Multiple R-squared:  0.1331, Adjusted R-squared:  0.1212 
## F-statistic: 11.14 on 2 and 145 DF,  p-value: 3.169e-05

# TRANSFORMACION SERIES AL DOMINIO DE LA FRECUENCIA
K <- c(rep(1,length(ipi)))
ipi.1 = gdf(ipi)
t.1=gdf(t)
t.2=gdf(t2)
K.1=gdf(K)

# REGRESIÓN EN EL DOMININO DE LA FRECUENCIA
summary(lm(ipi.1~0+K.1+t.1+t.2))

## 
## Call:
## lm(formula = ipi.1 ~ 0 + K.1 + t.1 + t.2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.195  -5.415  -0.469   2.552  65.236 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## K.1 95.2675702  2.9250060  32.570  < 2e-16 ***
## t.1  0.3413570  0.0906335   3.766  0.00024 ***
## t.2 -0.0025658  0.0005892  -4.355  2.5e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.7 on 145 degrees of freedom
## Multiple R-squared:  0.9872, Adjusted R-squared:  0.987 
## F-statistic:  3739 on 3 and 145 DF,  p-value: < 2.2e-16

Definiendo una matriz \(A\) de tamaño n x n de ceros excepto en las posiciones de la diagonal principal, correspondientes a las frecuencias que queremos incluir en la regresión, y premultiplicando \(\dot y\) y \(\dot X\) por \(A\) eleminamos determindas observaciones y las reemplazamos por ceros para realizar la regresión band spectrum. Devolver al dominio del tiempo a estas observaciones requiere:

\[y^* = W^{T}A\dot y = W^{T}AWy\] \[x^* = W^{T}A\dot x = W^{T}AWx \ \ \ \ \ \ (5)\]

Al regresar \(y^*\) sobre \(x^*\) obtenemos un \(\beta\) idéntico al estimador que obtendríamos al estimar por MCO \(\dot y\) frente a \(\dot x\).

Cuando se realiza la regresión band spectrum de esta manera, ocurre un problema asociado a los grados de libertad de la regresión de \(y^*\) sobre \(x^*\) que asumen los programas estadisticos convencionales, \(n - k\). Los grados de libertad reales serían únicamente \(n'- k\), donde \(n'\) es el numero de frecuencias incluidas en la regresión band spectrum.

Tan H.B and Ashley R (1999), señalan que el procedimiento de elaboración de una RBS consta de tres etapas:

1.- Transformar los datos originales del dominio del tiempo al dominio de la frecuencia utilizando series finitas de senos y cosenos. Implicaría premultiplicar los datos originales por la matriz ortogonal \(W\) sugerida por Harvey (1978).

2.- Permitir la variación a través de m bandas de frecuencia usando variables Dummy \(D_j^1,..D_j^m\). Estas variables se elaboran a partir de submuestras de las \(n\) observaciones del dominio de frecuencias. De esta forma, \(D_j^s=\dot x_{j,k}\) si la observación \(j\) está en la banda de frecuencias \(s\) y \(D_j^s=0\), en el resto de los casos.

3.- Re-estimar el resultado del modelo de regresión en el dominio del tiempo con las estimaciones y los coeficientes de las \(m\) variables Dummy. Implicaría premultiplicar la ecuación de regresión ampliada por las variables Dummy por la transpuesta de \(W\).

Asumiendo entonces que las series, \(y_t\), \(x_t\), \(\beta_t\) y \(u_t\) pueden ser transformadas en el dominio de la frecuencia:

\[y_t=\eta^y+\sum_{j=1}^N[a^y_j\cos(\omega_j)+b^y_j\sin(\omega_j)\]

\[x_t=\eta^x+\sum_{j=1}^N[a^y_j\cos(\omega_j)+b^y_j\sin(\omega_j)]\]

\[\beta_t=\eta^\beta+\sum_{j=1}^N[a^\beta_j\cos(\omega_j)+b^\beta_j\sin(\omega_j)]\]

\[(6)\]

Premultiplicando (6) por \(W\) obtenemos:

\[\dot y=\dot x\dot\beta+\dot u \ \ \ \ \ \ (7)\]

donde \(\dot y = Wy\),\(\dot x = Wx\), \(\dot \beta = W\beta\) y \(\dot u = Wu\)

El sistema (7) puede reescribirse como:

\[\begin{equation} \dot y=Wx_tI_nW^T\dot \beta + WI_nW^T\dot u \end{equation}\]

Si denominamos \(\dot e=WI_nW^T\dot u\), podrían buscarse los \(\dot \beta\) que minimizaran la suma cuadrática de los errores \(e_t=W^T\dot e\).

Una vez encontrada la solución a dicha optimización, se transformarían las series al dominio del tiempo para obtener el sistema (6).

Para obtener una solución a la minimización de los errores que ofrezca el mismo resultado que la regresión lineal por mínimos cuadrados ordinarios, requiere utilizar una matriz de regresores cuya primera columna sería el vector de tamaño \(N\), \((1,0,0,...)\), la segunda columna sería la primera fila de la matriz \(WI_NX_tW^T\), y las columnas siguientes, corresponderían las a las frecuencias de senos o cósenos que queremos regresar.

Denominando a esta nueva matriz, de tamaño (Nxp), \(\dot X\), donde \(p=2+j\), siendo \(j\) las frecuencias de seno y coseno elegidas como explicativas, los coeficientes de la solución MCO serían:

\[\dot \beta = (\dot X' \dot X)^{-1}\dot X' \dot y\]

donde \(\beta_{0,1}\) sería el parámetro asociado a la constante, \(\beta_{1,1}\) el asociado a la pendiente, y \(\beta_{1,j}\) los asociados a las frecuencias de senos y cósenos elegidas.

La regresión en el dominio de la frecuencia para el IPI de Cantabria realizado con un filtrado de altas frecuencias, se muestra aqui:

# CREAMOS LA MATRIZ DE REGRESORES FILTRANDO LAS ALTAS FRECUENCIAS CON LA MATRIZ DE HARVEY
n=length(ipi)
M <- MW(length(ipi))
z <- lm(ipi~t+t2)$fitted 
cx <- M%*%diag(z)
cx <- cx%*%t(M)
id <- seq(1,n)
S1 <- data.frame(cx)
S2 <- S1[1:(2+(n/12)),]
X <- as.matrix(S2)
X <- t(X)



# REGRESION FILTRANDO LAS ALTAS FRECUENCIAS
rbs.fit=lm(ipi.1~0+X)
summary(rbs.fit)

## 
## Call:
## lm(formula = ipi.1 ~ 0 + X)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.204  -4.326  -0.237   1.651  65.231 
## 
## Coefficients:
##     Estimate Std. Error t value Pr(>|t|)    
## X1  12.16502    0.10704 113.646  < 2e-16 ***
## X2  -0.01543    0.10793  -0.143   0.8865    
## X3   0.16244    0.10615   1.530   0.1283    
## X4   0.06540    0.10725   0.610   0.5430    
## X5  -0.44072    0.10683  -4.125 6.46e-05 ***
## X6  -0.05344    0.10711  -0.499   0.6186    
## X7   0.26225    0.10697   2.452   0.0155 *  
## X8   0.19221    0.10706   1.795   0.0748 .  
## X9  -0.12187    0.10701  -1.139   0.2568    
## X10 -0.08563    0.10703  -0.800   0.4251    
## X11  0.03593    0.10702   0.336   0.7376    
## X12  0.05775    0.10698   0.540   0.5902    
## X13  0.05649    0.10685   0.529   0.5979    
## X14 -0.04020    0.10680  -0.376   0.7072    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.86 on 134 degrees of freedom
## Multiple R-squared:  0.9898, Adjusted R-squared:  0.9888 
## F-statistic: 932.6 on 14 and 134 DF,  p-value: < 2.2e-16

plot(ts(ipi,frequency = 12, start = c(2002, 1)),type="l",main="IPI.Cantabria",ylab="")
lines(ts(lm(ipi~t+t2)$fitted,frequency = 4, start = c(1980, 1)),type="l",col=2)
lines(ts(gdt(rbs.fit$fitted),frequency = 4, start = c(1980, 1)),type="l",col=3)
legend("top", ncol=3,c("ipi","Estimado LM","Estimado RBS"),cex=0.6,bty="n",fill=c(1,2,3))

Una manera de obtener pronosticos en el caso de la regresión RBS, es generar los correspondientes armónicos y utilizar los coeficientes de la regresión RBS para componer el pronostico.

Por pasos:

1. Convertimos al dominio de tiempo cada regresor (cada columna de la matriz X)

2. Se reliza la estimación en el dominio del tiempo la regresión MCO \(X_{1t}=\hat \beta_0 + \hat \beta_1 t + \hat \beta_1 t^2\) y su pronóstico

3. Se reliza la estimación en el dominio del tiempo la regresión MCO \(X_{2t}=\hat \beta_0 + \hat \beta_1 cos( \frac {2*\pi t} {T})\) y su pronóstico

4. Se reliza la estimación en el dominio del tiempo la regresión MCO \(X_{3t}=\hat \beta_0 + \hat \beta_1 sin( \frac {2*\pi t} {T})\) y su pronostico

5. Se reliza la estimación en el dominio del tiempo la regresión MCO \(X_{4t}=\hat \beta_0 + \hat \beta_1 cos( \frac {4*\pi t} {T})\) y su pronóstico

6. Se reliza la estimación en el dominio del tiempo la regresión MCO \(X_{5t}=\hat \beta_0 + \hat \beta_1 sin( \frac {4*\pi t} {T})\) y su pronostico

7. y sucesivamente ….

8. Se combinan las series pronosticadas con los coeficientes que resultan de la regresión en el dominio de la frecuencia.

A continuación se realiza un ejemplo para el IPI de Cantabria, utilizando los dos primeros armónicos.

# REGRESION FILTRANDO LAS ALTAS FRECUENCIAS
S2 <- S1[1:5,]
X <- as.matrix(S2)
X <- t(X)
rbs.fit.2=lm(ipi.1~0+X)
summary(rbs.fit)

## 
## Call:
## lm(formula = ipi.1 ~ 0 + X)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.204  -4.326  -0.237   1.651  65.231 
## 
## Coefficients:
##     Estimate Std. Error t value Pr(>|t|)    
## X1  12.16502    0.10704 113.646  < 2e-16 ***
## X2  -0.01543    0.10793  -0.143   0.8865    
## X3   0.16244    0.10615   1.530   0.1283    
## X4   0.06540    0.10725   0.610   0.5430    
## X5  -0.44072    0.10683  -4.125 6.46e-05 ***
## X6  -0.05344    0.10711  -0.499   0.6186    
## X7   0.26225    0.10697   2.452   0.0155 *  
## X8   0.19221    0.10706   1.795   0.0748 .  
## X9  -0.12187    0.10701  -1.139   0.2568    
## X10 -0.08563    0.10703  -0.800   0.4251    
## X11  0.03593    0.10702   0.336   0.7376    
## X12  0.05775    0.10698   0.540   0.5902    
## X13  0.05649    0.10685   0.529   0.5979    
## X14 -0.04020    0.10680  -0.376   0.7072    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.86 on 134 degrees of freedom
## Multiple R-squared:  0.9898, Adjusted R-squared:  0.9888 
## F-statistic: 932.6 on 14 and 134 DF,  p-value: < 2.2e-16

# Tendencia cuadrática
T=148
n=12
X.1=gdt(X[,1])
z.1=c(lm(X.1~t+t2)$fitted,predict(lm(X.1~t+t2),data.frame(t=seq(from=T+1, to=T+n),t2=seq(from=T+1, to=T+n)^2)))
plot(z.1,type="l", col="red")
lines(X.1, type="p")

# Armonico 1
X.2.1=gdt(X[,2])
z.2.1=c(lm(X.2.1~cos(2*pi*t/T))$fitted,predict(lm(X.2.1~cos(2*pi*t/T)),data.frame(t=seq(from=T+1, to=T+n))))
plot(z.2.1,type="l", col="red")
lines(X.2.1, type="p")

X.2.2=gdt(X[,3])
z.2.2=c(lm(X.2.2~sin(2*pi*t/T))$fitted,predict(lm(X.2.2~sin(2*pi*t/T)),data.frame(t=seq(from=T+1, to=T+n))))
plot(z.2.2,type="l", col="red")
lines(X.2.2, type="p")

# Armónico 2
X.3.1=gdt(X[,4])
z.3.1=c(lm(X.3.1~cos(4*pi*t/T))$fitted,predict(lm(X.3.1~cos(4*pi*t/T)),data.frame(t=seq(from=T+1, to=T+n))))
plot(z.3.1,type="l", col="red")
lines(X.3.1, type="p")

X.3.2=gdt(X[,5])
z.3.2=c(lm(X.3.2~sin(4*pi*t/T))$fitted,predict(lm(X.3.2~sin(4*pi*t/T)),data.frame(t=seq(from=T+1, to=T+n))))
plot(z.3.2,type="l", col="red")
lines(X.3.2, type="p")

# Pronostico serie de dos armónicos
fit.2=rbs.fit.2$coefficients[1]*z.1+rbs.fit.2$coefficients[2]*z.2.1+rbs.fit.2$coefficients[3]*z.2.2+rbs.fit.2$coefficients[4]*z.3.1+rbs.fit.2$coefficients[5]*z.3.2
plot(ipi)
lines(fit.2,col="red")

Para obtener pronosticos con mayor precision, vamos a crea una serie de funciones que lo que hacer es generar indices para matriz ortogonal, \(W\), con \(m\) adelantos.

MW modificada

En la función MW2, m desfasa la serie original

MW2 <- function(n,m) {
# Author: Francisco Parra Rodriguez 
# Some ideas from: Harvey, A.C. (1978), Linear Regression in the Frequency Domain, International Economic Review, 19, 507-512.
# http://econometria.wordpress.com/2013/08/21/estimation-of-time-varying-regression-coefficients/ 
uno <- as.numeric (m:(n+m-1))
A <- matrix(rep(sqrt(1/n),n), nrow=1)
if(n%%2==0){
for(i in 3:n-1){ 
if(i%%2==0) {
A1 <- matrix(sqrt(2/n)*cos(pi*(i)*(uno-1)/n), nrow=1)
A <- rbind(A,A1)}
 else {
A2 <- matrix(sqrt(2/n)*sin(pi*(i-1)*(uno-1)/n), nrow=1)
A <- rbind(A,A2)
}} 
AN <- matrix(sqrt(1/n)*(-1)^(uno+1), nrow=1)
A <- rbind(A,AN)
A
} else {
for(i in 3:n-1){ 
if(i%%2==0) {
A1 <- matrix(
sqrt(2/n)*cos(pi*(i)*(uno-1)/n), nrow=1)
A <- rbind(A,A1)}
 else {
A2 <- matrix(sqrt(2/n)*sin(pi*(i-1)*(uno-1)/n), nrow=1)
A <- rbind(A,A2)
}} 
AN <- matrix(
sqrt(2/n)*sin(pi*(n-1)*(uno-1)/n), nrow=1)
A <- rbind(A,AN)
}
}

función gdf transformada

En la función gdf2, m adelanta la serie original

gdf2 <- function(y,m) {
  a <- matrix(y,nrow=1)
  n <- length(y)
  A <- MW2(n,m)
  A%*%t(a)
}

función gdt transformada

En la función gdt2, adelanta la serie original

gdt2 <- function(y,m) {
  # Author: Francisco Parra Rodriguez 
  # http://econometria.wordpress.com/2013/08/21/estimation-of-time-varying-regression-coefficients/ 
  a <- matrix(y,nrow=1)
  n <- length(y)
  A <- MW2(n,m)
  t(A)%*%t(a)
}

Utilizando estas nuevas matrices, la serie del IPI regional de Cantabria se pronosticaria con las altas frecuencias de esta forma:

# CREAMOS LA MATRIZ DE REGRESORES FILTRANDO LAS ALTAS FRECUENCIAS CON LA MATRIZ DE HARVEY
n=length(ipi)
M2 <- MW2(length(ipi),12)
z <- predict(lm(ipi~t+t2),data.frame(t=seq(from=11, to=n+10),t2=seq(from=11, to=n+10)^2))
cx <- M2%*%diag(z)
cx <- cx%*%t(M2)
id <- seq(1,n)
S1 <- data.frame(cx)
S2 <- S1[1:abs(2+(n/12)),]
X <- as.matrix(S2)
X <- t(X)



# PRONOSTICO FILTRANDO LAS ALTAS FRECUENCIAS
coef=t(as.matrix(rbs.fit$coefficients))
rbs.fit.3=coef%*%t(X)
dim(rbs.fit.3)

## [1]   1 148

fit.3=gdt2(rbs.fit.3,12)
plot(window(ts(ipi,frequency = 12, start = c(2002, 1)),start = c(2003,1)),type="l",main="IPI.Cantabria",ylab="")
lines(ts(z,frequency = 12, start = c(2003, 1)),type="l",col=2)
lines(ts(fit.3,frequency = 12, start = c(2003, 1)),type="l",col=3)
legend("top", ncol=3,c("ipi","Estimado LM","Estimado RBS"),cex=0.6,bty="n",fill=c(1,2,3))

La regresión en el dominio de la frecuencia tambien puede utilizarse para modelos bivariados. En la libreria descomponer, hay una función “rdf” para realizar dicha regresión.

El algoritmo de calculo “rdf” se realiza en las siguentes fases:

1. Calcula el co-espectro de la serie “x” e “y”

Sea \(x\) un vector n x 1 el modelo transformado en el dominio de la frecuencia esta dado por: \(\hat x= Wx\).

Sea \(y\) un vector n x 1 el modelo transformado en el dominio de la frecuencia esta dado por: \(\hat y= Wy\)

Denominando \(p_j\) el ordinal del cross-periodograma de \(\hat x\) y \(\hat y\) en la frecuencia \(\lambda_j=2\pi j/n\), y \(\hat x_j\) el j-th elemento de \(\hat x\) y \(\hat y_j\) el j-th elemento de \(\hat y\), entonces

\[ \left\lbrace \begin{array}{ll} p_j=\hat x_{2j}\hat y_{2j}+\hat x_{2j+1}\hat y_{2j+1} & \forall j = 1,...\frac{n-1}{2}\\ p_j=\hat x_{2j}\hat y_{2j}& \forall j = \frac{n}{2}-1 \end{array} \right. \]

\[p_0=\hat x_{1}\hat y_{1}\]

2. Ordena el co-espectro en base al valor absoluto de \(|p_j|\) y genera un índice en base a ese orden para cada coeficiente de fourier.

3. Calcula la matriz \(Wx_tI_nW^T\) y la ordena en base al indice anterior.

4. Obtiene \(\dot e=WI_nW^T\dot u\), incluyendo el vector correspondiente al parámetro constante, \((1,0,...0)^n\), y calucula el modelo utilizando los dos primeros regresores de la matriz \(Wx_tI_nW^T\) reordenada y ampliada, calcula el modelo para los 4 primeros, para los 6 primeros, hasta completar los \(n\) regresores de la matriz.

5. Realiza el Test de Durbin (1967) a los modelos estimados, y selecciona aquellos en donde los \(e_t=W^T\dot e\) están dentro de las bandas elegidas.

Denominando \(p_j\) al ordinal del periodograma de \(\hat a\) en la frecuencia \(\lambda_j=2\pi j/n\), y \(\hat a_j\) el j-th elemento de \(\hat a\), entonces

\[\left\lbrace \begin{array}{ll} p_j=\hat a_{2j}^{2}+\hat a_{2j+1}^{2} & \forall j = 1,...\frac{n-1}{2}\\ p_j=\hat a_{2j}^{2}& \forall j = \frac{n}{2}-1 \end{array} \right.\]

\[p_0=\hat a_{1}^{2}\]

Entonces, el cuadrado de \(\hat a\) puede ser utilizado como un estimador consistente del periodograma de \(a\).

El test de Durbin está basado en el siguiente estadistico:

\[s_j=\frac{\sum_{r=1} ^j p_r}{\sum_{r=1}^m p_r}\]

donde \(m=\frac{1}{2}n\) para \(n\) par y \(\frac{1}{2}(n-1)\) para \(n\) impar.

El estadístico \(s_j\) ha de encontrarse entre unos límites inferior y superior de valores críticos que han sido tabulados por Durbin (1969). Si bien, hay que tener presente que el valor \(p_0\) no se considera en el cálculo del estadístico, esto es, \(p_0=\hat v_1=0\).

6. De todos ellos, elige aquél que tiene menos regresores. Si no encuentra modelo, devuelve la solución MCO.

A continuación va a realizarse una RBS con datos trimestrales de la productividad y salarios reales en Canada, correspondientes al primer trimestre de 1980 finalizando el cuarto trimestre de 2000. En la libreria “descomponer”, las funciones “gdf” y “gdt” transforman las series del tiempo a la frecuencia y de la frecuencia al tiempo, siguiendo la transformación sugerida por Harvey (1978).

# LEEMOS DATOS
library(vars)
data("Canada")
head(Canada)

##                e     prod       rw    U
## 1980 Q1 929.6105 405.3665 386.1361 7.53
## 1980 Q2 929.8040 404.6398 388.1358 7.70
## 1980 Q3 930.3184 403.8149 390.5401 7.47
## 1980 Q4 931.4277 404.2158 393.9638 7.27
## 1981 Q1 932.6620 405.0467 396.7647 7.37
## 1981 Q2 933.5509 404.4167 400.0217 7.13

P <- as.numeric(Canada[, "prod"])
W <- as.numeric(Canada[, "rw"])

# REGRESIÓN MINIMO CUADRADA
summary(lm(W~P))

## 
## Call:
## lm(formula = W ~ P)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.485  -8.761  -1.274   7.755  29.873 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1242.4281   165.0193  -7.529 5.94e-11 ***
## P               4.1273     0.4046  10.200 3.00e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.54 on 82 degrees of freedom
## Multiple R-squared:  0.5593, Adjusted R-squared:  0.5539 
## F-statistic:   104 on 1 and 82 DF,  p-value: 3.003e-16

# TRANSFORMACION SERIES AL DOMINIO DE LA FRECUENCIA
K <- c(rep(1,84))
P.1 = gdf(P)
W.1=gdf(W)
K.1=gdf(K)

# REGRESIÓN EN EL DOMININO DE LA FRECUENCIA
summary(lm(W.1~0+K.1+P.1))

## 
## Call:
## lm(formula = W.1 ~ 0 + K.1 + P.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -86.668  -4.440  -3.030  -1.072  29.500 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## K.1 -1242.4281   165.0193  -7.529 5.94e-11 ***
## P.1     4.1273     0.4046  10.200 3.00e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.54 on 82 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 3.383e+04 on 2 and 82 DF,  p-value: < 2.2e-16

# REGRESION EN EL DOMINIO DE LA FRECUENCIA FILTRANDO LAS FRECUENCIAS RELEVANTES 
rdf.mod=rdf(W,P)
str(rdf.mod)

## List of 6
##  $ datos      :'data.frame': 84 obs. of  4 variables:
##   ..$ Y  : num [1:84] 386 388 391 394 397 ...
##   ..$ X  : num [1:84] 405 405 404 404 405 ...
##   ..$ F  : num [1:84] 395 391 387 391 401 ...
##   ..$ res: num [1:84] -9.31 -2.58 3.85 2.57 -3.98 ...
##  $ Fregresores: num [1:84, 1:16] 1 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:84] "X1" "X2" "X3" "X4" ...
##   .. ..$ : chr [1:16] "C" "1" "2" "3" ...
##  $ Tregresores: num [1:84, 1:16] 0.109 0.109 0.109 0.109 0.109 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:16] "C" "1" "2" "3" ...
##  $ Nregresores: num 16
##  $ sse        : num 764
##  $ gcv        : num 13.9

plot(ts(W,frequency = 4, start = c(1980, 1)),type="l",main="Salarios reales.Canada",ylab="")
lines(ts(lm(W~P)$fitted, frequency = 4, start = c(1980, 1)), type="l", col=2)
lines(ts(rdf.mod$datos$F, frequency = 4, start = c(1980, 2)), type="l", col=3)
legend("top", ncol=3, c("W","Estimado LM","Estimado RBS"), cex=0.6, bty="n", fill=c(1,2,3))

Árbol de regresión y random forest.

La mayoría de los métodos de ML, no tienen conciencia del tiempo. Por el contrario, consideran que las observaciones son independientes e idénticamente distribuidas. Obviamente, esta suposición se viola en los datos de series de tiempo que se caracterizan por la dependencia en serie.

Los métodos basados en el árbol de decisión, por ejemplo, no pueden predecir una tendencia, es decir, no se extrapolan. Los árboles operan mediante reglas que dividen recursivamente el espacio de entrada. Por lo tanto, no pueden predecir valores que se encuentren fuera del rango de valores del objetivo en el conjunto de entrenamiento. Pero utilizando transformaciones y diferenciaciones pueden hacerse pronósticos de series de tiempo con bosques aleatorios (https://www.statworx.com/de/blog/time-series-forecasting-with-random-forest/).

Esencialmente, una serie de tiempo (univariante) es un vector de valores indexados por tiempo. Para que sea “aprendible”, necesitamos hacer un preprocesamiento. Esto puede incluir algunos o todos los siguientes:

• Transformaciones estadísticas (transformación Box-Cox, transformación logarítmica, etc.) • Detrending (diferenciación, STL, etc.) • Funciones de retardo de tiempo • Ingeniería de características (retrasos, estadísticas continuas, términos de Fourier, variables de tiempo, etc.)

Vamos a utilizar la serie de pasajeros en lineas aereas para realizar este ejercicio.. El pronostico de una serie temporal mediante estas técnicas requiere transformar la serie hasta hacerla estacionaria. Con nsdiffs elegimos el orden de diferenciación para hacerla estacionaria en media, y tomamos el logaritmo para que sea estacionaria en varianza.

ipi<-ts(ipi,c(2000,1),frequency=12)
x<-window(ipi, end = c(2010, 12))


# Buscamos el mejor orden de diferenciación
n_diffs <- nsdiffs(x)

# Realizamos una transformación logaritmica
x <- x %>% 
  log() %>% 
  diff(n_diffs)

# Respresentamos los resultados
plot(x)

Construimos una matriz con la función R “embed”. Cada fila de la matriz resultante consta de secuencias x [t], x [t-1], …, x [t-dimension + 1], donde t es el índice original de x. Si x es una matriz, es decir, x contiene más de una variable, entonces x [t] consiste en la tth observación en cada variable. En en ejemplo la primera columna empieza con el dato de Julio de 2000 (seis meses adelante), la siguiente con el dato de junio de 2000, y así sucesivamente. En consecuencia, en la ultima estaría la serie original.

lag_order <- 6 # numero de meses de desfase (seis meses)
horizon <- 12 # the forecast horizon (12 meses)

x_mbd <- embed(x, lag_order + 1) # embedding magic!

Seleccionamos dos muestras, en la de entrenamiento están los de 2000 a 2010, y en la de test los de 2011.

y_train <- x_mbd[, 1] # the target
X_train <- x_mbd[, -1] # everything but the target

y_test <- window(ipi, start = c(2011, 1), end = c(2011, 12)) # the year 2012
X_test <- x_mbd[nrow(x_mbd), c(1:lag_order)] # the test set consisting
# of the six most recent values (we have six lags) of the training set. It's the
# same for all models.

Estimamos un modelo con randomForest.

library(randomForest)

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

forecasts_rf <- numeric(horizon)

for (i in 1:horizon){
  # set seed
  set.seed(2019)

  # Estimo el modelo
  fit_rf <- randomForest(X_train, y_train)

  # Realizo los prnosticos para la muestra test
  forecasts_rf[i] <- predict(fit_rf, X_test)

  
# Reestructuramos repetidamente los datos de entrenamiento para reflejar la distancia de tiempo correspondiente al horizonte de pronóstico actual.
  y_train <- y_train[-1] 

  X_train <- X_train[-nrow(X_train), ] 
}

Dada la transformación logarítmica realizada, se realiza la transformación inversa. Es decir, primero revertimos la diferenciación y luego la transformación logarítmica. Hacemos esto exponiendo la suma acumulativa de nuestros pronósticos transformados y multiplicando el resultado con la última observación de nuestra serie de tiempo. En otras palabras, calculamos:

$ y_{t+H} = y_{t}exp(\sum_{h = 1}^H \tilde{y}_{t+h}\right)$

# exponenciamos la suma acumulativa
exp_term <- exp(cumsum(forecasts_rf))

# extraemos la ultima obsrvación de nuestra base de datos (y_t)
last_observation <- as.vector(ipi[144])

# calculamos las prdiciones resultantes
backtransformed_forecasts <- last_observation * exp_term

# convertimos a formato ts y lo representamos
y_pred <- ts(
  backtransformed_forecasts,
  start = c(2012, 1),
  frequency = 12
)

IPI.f=ts(c(ipi,y_pred),frequency = 12,start=c(2000,1)) 
  
# Representamos la serrie con los resultados estimados
plot(IPI.f)

Redes Neuronales Artificiales

La RNA es una técnica de inferencia basada en Machine Learnig (aprendizaje máquina). Existen distintos modelos de redes neuronales, los más utilizados son el “Perceptron” y el “backpropagation”.

Hay que tener presente que una red neuronal se considera una “caja negra”, donde lo importante es la predicción, y no cómo se hace.

El proceso incluye, como ocurre con este tipo de técnicas, una fase de entrenamiento (training) para la optimización de las predicciones, y otra de test o prueba.

Los elementos de la red son:

• Las neuronas o nodos
• Las capas
  • De entrada
  • De salida
  • Oculta (puede tener a su vez varias capas)
• Los pesos
• La función de combinación
• La función de activación
• El objetivo (target)

Los nodos (neuronas) de la capa de entrada, se combinan con los nodos de la capa oculta mediante la función de combinación, que suele ser una combinación lineal de los nodos de entrada mediante los pesos. A las neuronas de las capas ocultas, se les aplica una función de activación, que suele ser la tangente hiperbólica de la anterior combinación más un parámetro por nodo oculto, con la que estimamos las neuronas de la capa de salida, y sus errores.

El siguiente script utiliza el package neuralnet para hacer una regresión, utilizando el ejemplo del IPI utilizados en el ejemplo anterior. En lo relativo al modelo de red neuronal que se va a estimar, destacar:

• El parametro threshol = 0.05 indica que las iteraciones se detendrán cuando el “Cambio” del error sea menor a 5% entre una iteracion de optimizacion y otra. Este “Cambio” es calculado como la derivada parcial de la funcion de error respecto a los pesos.

• El parámetro algorithm = “rprop+” refiere al algoritmo “Resilient Backpropagation”, que actualiza los pesos considerando únicamente el signo del cambio, es decir, si el cambio del error es en aumento (+) o disminución (-) entre una iteración y otra. Para detalles ver: https://en.wikipedia.org/wiki/Rprop

• El parámetro hidden = 1 especifica una capa oculta con 1 neuronas.

#LIBRERIAS Y DATOS
library(MASS); library(neuralnet); library(ggplot2)

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:randomForest':
## 
##     margin

set.seed(65)

 
# NORMALIZACION DE VARIABLES
datos_nrm <- as.data.frame(x_mbd)
names(datos_nrm)=c("x1","x2","x3","x4","x5","x6","IPI")
train_nrm <- datos_nrm[1:113,  ]
test_nrm  <- datos_nrm[113:125,, ]
 
# FORMULA
nms  <- c("x1","x2","x3","x4","x5","x6","IPI")
frml <- as.formula(paste("IPI ~", paste(nms[!nms %in% "IPI"], collapse = " + ")))
 
# MODELO
modelo.nn <- neuralnet(frml,data=train_nrm,threshold = 0.05,hidden=1,algorithm = "rprop+")                                             
 
# PREDICCION
pr.nn   <- compute(modelo.nn,within(test_nrm,rm(IPI)))
 
# Red
plot(modelo.nn)

# exponenciamos la suma acumulativa
exp_term <- exp(cumsum(pr.nn$net.result))

# extraemos la ultima obsrvación de nuestra base de datos (y_t)
last_observation <- as.vector(ipi[144])

# calculamos las prdiciones resultantes
backtransformed_forecasts <- last_observation * exp_term

# convertimos a formato ts y lo representamos
y_pred <- ts(
  backtransformed_forecasts,
  start = c(2012, 1),
  frequency = 12
)

IPI.f=ts(c(ipi,y_pred),frequency = 12,start=c(2000,1)) 
  
# Representamos la serrie con los resultados estimados
plot(IPI.f)

Agrupamiento de series temporales con DTW

Las técnicas de machine leaning tambien pueden aplicarse a la clasificación o agrupación de las series temporales en base a sus características comunes. Se trataría de dividir el data-set de series temporales en grupos basados en similitudes o distancias, de modo que las series temporales de un mismo grupo tengan la misma naturaleza dinámica. Para la agrupación de series de tiempo, el primer paso es elaborar una métrica de distancia ó similitud apropiada, y luego, en el segundo paso, utilizar técnicas de agrupación existentes, como k-means, clústering jerárquico, clústering basado en densidad o clústering subespacial.

El objetivo es formar grupos homogéneos de series temporales, con una similitud mínima entre grupos y máxima entre grupos. Para lo cual precisamos:

1. Medida de distancia

2. Prototipo (resume las características de todas las series en un grupo)

3. Algoritmo de clúster (particional o jerárquico más común)

Para encontrar estructuras de agrupamiento, de series temporales la metrica más utilizada es Dynamic Time Warping (DTW) (Sakoe & Chiba, 1978). Esta distancia es una técnica desarrollada en el ambito de reconocimiento de la voz, cuyo objetivo es comparar dos emisiones vocales aisladas a fin de determinar si corresponden o no a la misma palabra. El DTW se diseño para comprobar si existe una sincronización temporal (alineamiento temporal) entre dos grupos fónicos. La falta de alineamiento, por lo general, no obedece a una ley fija (p. e., un retardo constante), sino que se da de forma heterogénea, produciéndose así variaciones localizadas que aumentan o disminuyen la duración del tramo de análisis.

El problema de calcular una función de alineamiento se resuelve en el DTW mediante técnicas de programación dinámica, motivo por el que se conoce a esta técnica como Alineamiento Temporal Dinámico.

En general, el funcionamiento del algoritmo se basa en la búsqueda de un camino óptimo de coincidencia entre dos secuencias con ciertas restricciones. Las secuencias son transformadas no linealmente en el tiempo, de modo que son comprimidas o expandidas en el tiempo para que tengan el mismo largo, y así poder compararlas punto a punto.

De forma general, el algoritmo DTW consiste en una serie de pasos que se explican a continuación:  Se requieren dos matrices para el cálculo:

D  matriz de distancias.

G  matriz de programación dinámica, para obtener el camino mínimo entre el punto final y el inicial.

Se parte de 2 vectores: \(A= [a_1, a_2, a_3, .. , a_i , .., a_M]\) y \(B= [b_1, b_2, b_3, .. , b_j , .. , b_N]\)

Se denota como \(C= [c(1), c(2), .., c(k), .., c(K)]\) a la función de alineamiento, donde cada \(c(k)= [i(k) j(k)]\)tiene una matriz de coste \(d[c(k)]\), que muestra la diferencia entre los elementos comparados.

La función de coste tomando como medida la distancia euclieda, será \(d[c(k)]= (a_{i(k)}-b_{j(k)})^2\) y la función de alineamiento que minimice el coste total será:

\[D(C)=\sum_{k=1}^K d[c(k)] w(k)\]

donde \(W_k\) es un coeficiente de ponderación.

La función de alineamiento debe cumplir una serie de restricciones:

1. Monotonicidad: El camino debe tener el sentido de izquierda a derecha, y de arriba debajo de la matriz G.

2. Continuidad: El paso de un nodo al siguiente debe hacerse de modo que no haya nodos intermedios.

3. Frontera: El principio y el final deben coincidir con ciertos puntos, independientemente del camino seguido.

4. Ventana: Condición de ajuste (límite global): \(|i(k)-j(k)|\leq \Delta\), siendo \(\Delta\) un entero positivo.

5. Pendiente: Se aplican restricciones (límites locales) sobre las pendientes del camino para encontrar el camino óptimo.

y una serie de restricciones locales ((Sakoe & Chiba,1978).

El algoritmo DTW puede construirse con una medida de distancia eucliedea u otra medidad de distancia (por ejemplo,manhatan)

El principal problema de esta medida de similitud es el gran coste computacional que supone, ya que al tratarse de programación dinámica dadas dos cadenas de longitud \(N\) y \(M\) respectivamente el coste computacional del algoritmo es del orden de \(N x M\).

En el ejemplo, se utiliza un conjunto de datos de la serie temporal de gráficos de control sintético, que contiene 600 ejemplos de gráficos de control. Cada gráfico de control es una serie de tiempo con 60 valores. Hay seis clases:

1. 1-100 Normal,
2. 101-200 Cíclico,
3. 201-300 Tendencia creciente,
4. 301-400 Tendencia decreciente,
5. 401-500 Cambio hacia arriba y
6. 501-600 Descenso cambio.

Los datos y el script se ha obtenido de :http://www.rdatamining.com/examples/time-series-clústering-classification.

# lectura de datos
sc <- read.table("http://kdd.ics.uci.edu/databases/synthetic_control/synthetic_control.data", header=F, sep="")
# Selección aleatoria de variables
n <- 10
s <- sample(1:100, n)
idx <- c(s, 100+s, 200+s, 300+s, 400+s, 500+s)
sample2 <- sc[idx,]
observedLabels <- c(rep(1,n), rep(2,n), rep(3,n), rep(4,n), rep(5,n), rep(6,n))
# calculo de distancias DTW
library(dtw)

## Loading required package: proxy

## 
## Attaching package: 'proxy'

## The following object is masked from 'package:Matrix':
## 
##     as.matrix

## The following objects are masked from 'package:stats':
## 
##     as.dist, dist

## The following object is masked from 'package:base':
## 
##     as.matrix

## Loaded dtw v1.20-1. See ?dtw for help, citation("dtw") for use in publication.

distMatrix <- dist(sample2, method="DTW")
# clúster jerrárquico
hc <- hclust(distMatrix, method="average")
plot(hc, labels=observedLabels, main="",cex=0.5)