Estadística y Machine Learning con R

Inferencia no parámetrica

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, pueden 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.

\[\begin{equation} y_i= \beta_1 + \beta_2 x_i + e_i, i=1,2,…, n \end{equation}\]

.

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

Este es un modelo estadístico con tres parámetros desconocidos:\(\beta_1\); \(\beta_2\) 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_1, \beta_2)\) y \(m(x_i;\theta)=\beta_1 + \beta_2 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 de 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 nucleo).

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 una 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 o amplitudes de banda (windows, bandwidths) que deben tender a cero ”lentamente" (\(h \Rightarrow 0\) ,\(nh\Rightarrow inf\) ) 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

\(\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\) asignado los pesos obtenidos mediante la función núcleo se resuelve 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 función de regresión. Valores pequeños de 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 ajustes a los datos pero predicican mejor, pero si h es demasiado elevado tendremos una falta de ajuste a los datos (sub-estimación).

Utilizando la base de datos “cars” de R, que contine las variables “dist” (distancia de parada) y “speed” (velocidad), vamos a realizar la representación gráfica de la regresión kernel realizada con el estimador de Nadaraya–Watson con diferentes parámetros de suavizado.

data(cars)
plot(cars$speed, cars$dist)
lines(ksmooth(cars$speed, cars$dist, "normal", bandwidth = 2), col = 2)
lines(ksmooth(cars$speed, cars$dist, "normal", bandwidth = 5), col = 3)

Splines

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 bases de funciones para aproximar curvas con formas complicadas.

Las bases de spilines más populares:

• Bases de polinomios truncados.

• Bases de splines cúbicos.

• Bases de B-splines.

• Bases de thin plate splines.

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, 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 ]\).

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.

Splines de grado 0

Los splines de grado \(0\), se define en un solo tramo de nudo y ni siquiera es continua en los nudos. 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:

Splines lineal

Las funciones de spilines más comúnmente utilizadas son las de grado 3 ó 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 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_0 & 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.

Una de las bases de splines cúbicos más utilizadas basadas en \(q-2\) nodos interiores, \(x^*_j\), \(j=1,2,...,q-2\) es:

\[\begin{equation} \left\lbrace\begin{array}{ll} s_0(x)=1 \\ s_1(x)=x \\ s_{j+2}(x)=R[x_j,x^*_j] \end{array}\right.\end{equation}\]

Siendo \(R[x,z]=\frac 1 4[(z-\frac 1 2)^2-\frac 1 {12}][(x-\frac 1 2)^2-\frac 1 {12}]-\frac 1 {24}[(\mid z- x \mid - \frac 1 2)^4- [(\mid z- x \mid - \frac 1 2)^2+\frac 7 {240}]\)

Con esta base de splines definimos \(f\) a través de un modelo lineal con matriz de regresores \(X\) con \(n\) filas y \(q\) columnas cuya i_esima fila es:

\(X_i=[1,x_i,R[x_i,x^*_1],R[x_i,x^*_2],...,R[x_i,x^*_{q-2}]]\)

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

Splines lineal

Un tema importante es la elección del grado de suavización del spline. Una de las posibilidades es a través del contraste de hipótesis, valorar la posibilidad de utilizar uno o más nodos. Pero lo aconsejado es mantener fija la base de splines y controlar el grado de suavización añadiendo una penalización a la función objetivo de mínimos cuadrados:

\(\lambda \beta' S \beta\)

Donde \(S\) es una matriz de orden \(q*q\) con coeficientes conocidos que dependen de la base elegida y un parámetro de suavizado \(\lambda\).

La solución del modelo de regresión lineal penalizado en donde la matriz de regresores está ahora definida por la base de splines y la penalización sería:

\(\hat \beta = (X'X - \lambda S)^{-1}X'y\)

El modelo de regresión lineal con spilines penalizados es equivalente al siguiente modelo de regresión lineal:

\(Y=\beta X + e\)

En donde \(Y'=(y,0,0...0)\) es un vector de dimensión \((n+q)\), es decir el vector \(y\) seguido de tantos ceros como nodos se han utilizado en la base de los splines.

La matriz de regresores

\[X' = \begin{bmatrix}X\\ B \sqrt \lambda\\ \end{bmatrix}\]

tiene ahora orden \((n+q)*q\), siendo \(B\) una matriz que cumple \(B=S'S\) y que se obtiene a través de la descomposición de Cholesky, \(\lambda\) el parámetro de suavizado y \(e\) un vector de \((n+q)\) errores aleatorios.

El parámetro de suavización,\(\lambda\),es a priori desconocido y hay que determinarlo, si es muy alto suaviza los datos en exceso, un criterio utilizado para elegir el parámetro es del valor que minimiza el estadístico general de validación cruzada.

La regresión por splines puede realizarse con múltiples variables explicativas, si tenemos ahora dos explicativas, \(x_i\) y \(z_i\),y queremos estimar el siguiente modelo aditivo:

\(y_i=f_1(x_i)+f_2(z_i)+e_i\)

Representaríamos cada una de estas dos funciones a través de una base de splines penalizados, que tomando la base cúbica quedaría:

\(f_1(x)=\delta_1+\delta_2X_i+\sum_{j=1}^{q-2} \delta_jR[x_i,x_j^*]\)

\(f_2(z)=\gamma_1+\gamma_2X_i+\sum_{j=1}^{q-2} \gamma_jR[z_i,z_j^*]\)

Partiendo de la base de datos “cars”, la función R “smooth.spline” realiza la regresión por splines utilizando una base de splinee cúbicos penalizados:

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")

En la función “smooth.spline” el parámetro de suavizado es un valor generalmente entre 0 y 1, en tanto que el coeficiente \(\lambda\) se obtiene en el criterio de aceptación (logaritmo de verosimilitud penalizado). En el ejercicio el programa elige un \(\lambda=0.7801305\). Si se desea un función menos suavizada habrá que elegir un parámetro de suavizado más bajo, en línea roja se representa en el gráfico la regresión por splines que se obtendría con un parámetro de suavizado de valor \(0.10\).

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). \(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 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: splines

## Loading required package: Matrix

## 
## Attaching package: 'fda'

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

library(fda.usc)

## Loading required package: MASS

## Loading required package: mgcv

## Loading required package: nlme

## This is mgcv 1.8-17. 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")

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)

La función: smooth.basis(argvals=1:n, y, fdParobj), del R-package fda, donde 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")

Regresion 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 pasaser 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[a_j\cos(2\pi\frac{ft}n)+b_j\sin(2\pi\frac{ft}n)]\) (1)

donde \(\eta\) es la media de la serie, \(a_j\) y \(b_j\) son su amplitud,\(f\) son las frecuencias que 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). Las 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, asi :

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

El analisis espectral puede utilizarse para identificar y cuantificar en procesos aparentemente aperiodicos, sucesiones de cicos de periodo 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 carcter cíclico la hacer aparecer como un serie no periodica e irregular. De hecho la ecuación (2), muestra que cada observación \(t\) de una serie de tiempo, es el resultado 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 termino de error.

Una manera practica de pasar desde el dominio del tiempo al dominio de la frecuencia es pre-multiplicando 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).

\[\begin{equation} y=X\beta+u \end{equation}\] (4)

donde \(X\) es una matriz \(n x k\) de observaciones de \(k\) variables independientes, \(\beta\) es un vecto \(k x I\) de parámetros, \(y\) es un vecto \(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 autocorrelación en los residuos. Engle (1974) muesta que si los residuos están correlacionados serialmente y son generados por un procieso estacionario estocastico, 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 (4) 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 (4) por \(W\), obtenemos: \[\begin{equation} \dot y=\dot X\beta+\dot u \end{equation}\] (5)

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

Si el vector de las perturbaciones en (4) cumple las hipoteis 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 (5) 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 identico al estimador MCO de (4).

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

Definiendo una matriz de tamaño \(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 estas observaciones requiere:

\[\begin{equation} y^* = W^{T}A\dot y = W^{T}AWy \\ x^* = W^{T}A\dot x = W^{T}AWx \end{equation}\] (6)

Al regresar \(y^*\) sobre \(x^*\) obtenemos un \(\beta\) identico 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\), en vez de los grados de libertad reales que serían unicamente \(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 una matriz ortogonal, \(W\), sugerida por Harvey (1978).

2.- Permitir la variación de 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\).

Realizamos una RBS con datos trimestrales del PIB y empleo en Canada, correspondientes al primer trimestre de 1980 al cuarto trimestre del 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).

library(descomponer)

## Warning: package 'descomponer' was built under R version 3.4.1

## Loading required package: taRifx

library(vars)

## Loading required package: strucchange

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: urca

## Loading required package: lmtest

data("Canada")
summary(Canada)

##        e              prod             rw              U         
##  Min.   :928.6   Min.   :401.3   Min.   :386.1   Min.   : 6.700  
##  1st Qu.:935.4   1st Qu.:404.8   1st Qu.:423.9   1st Qu.: 7.782  
##  Median :946.0   Median :406.5   Median :444.4   Median : 9.450  
##  Mean   :944.3   Mean   :407.8   Mean   :440.8   Mean   : 9.321  
##  3rd Qu.:950.0   3rd Qu.:410.7   3rd Qu.:461.1   3rd Qu.:10.607  
##  Max.   :961.8   Max.   :418.0   Max.   :470.0   Max.   :12.770

E <- as.numeric(Canada[, "e"])
PIB <- as.numeric(Canada[, "prod"])-E
summary(lm(E~PIB))

## 
## Call:
## lm(formula = E ~ PIB)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.4951 -2.3631  0.5435  2.0502  7.1228 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 158.77705   32.62608   4.867 5.42e-06 ***
## PIB          -1.46426    0.06082 -24.077  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.243 on 82 degrees of freedom
## Multiple R-squared:  0.8761, Adjusted R-squared:  0.8746 
## F-statistic: 579.7 on 1 and 82 DF,  p-value: < 2.2e-16

K <- c(rep(1,84))
PIB.1 = gdf(PIB)
E.1=gdf(E)
K.1=gdf(K)
summary(lm(E.1~0+K.1+PIB.1))

## 
## Call:
## lm(formula = E.1 ~ 0 + K.1 + PIB.1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.5629  -0.7838  -0.1334   0.4203  18.2778 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## K.1   158.77705   32.62608   4.867 5.42e-06 ***
## PIB.1  -1.46426    0.06082 -24.077  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.243 on 82 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 3.56e+06 on 2 and 82 DF,  p-value: < 2.2e-16

Teniendo presente el periodograma del PIB y el empleo:

gperiodograma(PIB)

gperiodograma(E)

Utilizamos variables Dummys para separar los ciclos de menor frecuencia del resto, y realizamos una RBS para los ciclos de menor frecuencia.

D1 <- c(rep(1,10),rep(0,74)) 
PIB.1.1=D1*PIB.1
rbs.fit=lm(E.1~0+K.1+PIB.1.1)
plot(ts(E,frequency = 4, start = c(1980, 1)),type="l")
lines(ts(lm(E~PIB)$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("E","Estimado LM","Estimado RBS"),cex=0.6,bty="n",fill=c(1,2,3))

La regresión RBS cuando se utiliza diferenciando o filtrando frecuencias mediante variables Dummys, opera como un filtro frecuencial. Hay que tener presente que las series originales se transformas del tiempo a la frecuencia por la matriz ortogonal que hemos denominado \(W\), y dicha matriz tiene una dimensión (n X n), esto determina que la transformación de los datos de la matriz de regresores para la predicción, ahora de tamaño (n+1 x k), de ser transformados al dominio de la frecuencia por una matriz \(W\) de tamaño (n+1 x n+1), se encontrarían en una frecuencia distinta a la que precisa la separación que ser realiza con las variables Dummys, para solucionar estos problemas de falta de sintonía entre las series de frecuiencia, la predicción ha de realizarse sobre la matriz \(W\) de dimensión (n x n), y las observaciones (2 a n+1) de la matriz de regresores.

Vamos a realizar la predicción correspondiente a un incremento trimestral del empleo del 0,25%. Como la preducción se realiza en el dominio de la frecuencia hay transformar los datos al dominio del tiempo

PIB_2.1980_1.2001=c(PIB[2:84],PIB[84]*1.0025)
D1 <- c(rep(1,10),rep(0,74)) 
PIB.2 = gdf(PIB_2.1980_1.2001)
PIB.2.1=D1*PIB.2
E.2.fitted=rbs.fit$coefficients[1]*K.1+rbs.fit$coefficients[2]*PIB.2.1
E_2.1980_1.2001.fitted=gdt(E.2.fitted)
E_2.1980_1.2001.fitted[84]

## [1] 943.7846

plot(ts(E,frequency = 4, start = c(1980, 1)),type="l")
lines(ts(lm(E~PIB)$fitted,frequency = 4, start = c(1980, 1)),type="l",col=2)
lines(ts(E_2.1980_1.2001.fitted,frequency = 4, start = c(1980, 2)),type="l",col=3)
legend("top", ncol=3,c("E","Estimado LM","Estimado RBS"),cex=0.6,bty="n",fill=c(1,2,3))

Regresión en el dominio de la frecuencia con un modelo no lineal.

Asumiendo que las series, \(y_t\),\(x_t\),\(\beta_t\) and \(ut\), 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)]\]

Otenemos dichas series pre-multiplicando (7) por \(W\)

\(\dot y=\dot x\dot\beta+\dot u\) (8)

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

El sistema (8) 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 (7).

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 nueva esta 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.

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 a el 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 ampliadas, 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 () a los modelos estimados, y selecciona aquellos en donde los \(e_t=W^T\dot e\) están dentro de las bandas elegidas a los niveles de significación que determina la función: \(\alpha=0.1\).

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

Denominando \(p_j\) el 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 del \(\hat a\) puede ser utilizado como un estimador consistente del periodograma de \(a\).

El test de Durbin esta 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 en 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_o\) no se considera en el cálculo del estadístico esto es, \(p_o=\hat v_1=0\)

rdf.mod=rdf(E,PIB)
str(rdf.mod)

## List of 4
##  $ datos      :'data.frame': 84 obs. of  4 variables:
##   ..$ Y  : num [1:84] 930 930 930 931 933 ...
##   ..$ X  : num [1:84] -524 -525 -527 -527 -528 ...
##   ..$ F  : num [1:84] 930 931 932 932 931 ...
##   ..$ res: num [1:84] -0.251 -0.71 -1.344 -0.308 1.368 ...
##  $ Fregresores: num [1:84, 1:18] 1 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:84] "X1" "X2" "X3" "X4" ...
##   .. ..$ : chr [1:18] "C" "1" "2" "3" ...
##  $ Tregresores: num [1:84, 1:18] 0.109 0.109 0.109 0.109 0.109 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:18] "C" "1" "2" "3" ...
##  $ Nregresores: num 18

plot(ts(E,frequency = 4, start = c(1980, 1)),type="l")
lines(ts(lm(E~PIB)$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("E","Estimado LM","Estimado RBS"),cex=0.6,bty="n",fill=c(1,2,3))

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 el del “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 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 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.

Red Neuronal Artificial

El siguiente script utiliza el package neuralnet para hacer una regresion y predecir el valor de viviendas (expresado en $1000), usando el data set Boston incluido en el package MASS. Para un mejor ajuste del modelo (es decir, para que la red aprenda mejor), se hace un preprocesamiento de los datos, donde se normaliza el datas set Boston. En consecuencia, una vez se obtiene la predicción con la red neural, al presentarse en valores normalizados, hay que transformar los datos para obtener la valoración en términos de precios de la vivienda.

Como es habitual en los ejercicios de minería de datos, el conjunto de datos se divide en una muestra de entrenamiento (el 70%) y otra de test.

La programación procede de : http://apuntes-r.blogspot.com.es/2015/12/regresion-con-red-neuronal.html, en donde se destacan las siguientes notas:

• El parametro threshol = 0.05 indica que las iteraciones se detendran 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 parametro 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 parametro hidden = c(7,5) especifica una primera capa oculta con 7 neuronas y una segunda capa oculta con 5 neuronas.

#LIBRERIAS Y DATOS
# -----------------------------------------------------
library(MASS); library(neuralnet); library(ggplot2)
set.seed(65)
datos    <- Boston
n        <- nrow(datos)
muestra  <- sample(n,n*.70)
train    <- datos[muestra, ]
test     <- datos[-muestra, ]
 
 
# NORMALIZACION DE VARIABLES
# -----------------------------------------------------
maxs      <- apply(train, 2, max)
mins      <- apply(train, 2, min)
datos_nrm <- as.data.frame(scale(datos, center = mins, scale = maxs - mins))
train_nrm <- datos_nrm[muestra, ]
test_nrm  <- datos_nrm[-muestra, ]
 
 
# FORMULA
# -----------------------------------------------------
nms  <- names(train_nrm)
frml <- as.formula(paste("medv ~", paste(nms[!nms %in% "medv"], collapse = " + ")))
 
 
# MODELO
# -----------------------------------------------------
modelo.nn <- neuralnet(frml,
                       data          = train_nrm,
                       hidden        = c(7,5), # ver Notas para detalle 
                       threshold     = 0.05,   # ver Notas para detalle
                       algorithm     = "rprop+" 
                       )
 
 
# PREDICCION
# -----------------------------------------------------
pr.nn   <- compute(modelo.nn,within(test_nrm,rm(medv)))
 
# se transforma el valor escalar al valor nominal original
medv.predict <- pr.nn$net.result*(max(datos$medv)-min(datos$medv))+min(datos$medv)
medv.real    <- (test_nrm$medv)*(max(datos$medv)-min(datos$medv))+min(datos$medv)
 
 
 
# SUMA DE ERROR CUADRATICO
# -----------------------------------------------------
(se.nn <- sum((medv.real - medv.predict)^2)/nrow(test_nrm))

## [1] 8.730475173

#GRAFICOS
# -----------------------------------------------------
# Errores
qplot(x=medv.real, y=medv.predict, geom=c("point","smooth"), method="lm", 
      main=paste("Real Vs Prediccion. Summa de Error Cuadratico=", round(se.nn,2)))

## Warning: Ignoring unknown parameters: method

# Red
plot(modelo.nn)

Red Neuronal Artificial Precio Vivienda

Árbol de regresión

El Arbol de decisión puede ser utilizado para clasificar unidades y tambien para realizar una regresión no paramétricas (Árbol de regresión). Si utilizamos la función “rpart” hay que seleccionar en “method” la elección “ANOVA”. En la opción “control” encontramos parámetros opcionales para controlar el crecimiento de los árboles. Por ejemplo, en control = rpart.control (minsplit = 30, cp = 0.001) se le indica que el número mínimo de observaciones en un nodo sea 30 antes de intentar una división y que una división debe pararse cuando el Coste factor de complejidad sea inferior a un 0.001, con el método anova, esto significa que el R-cuadrado total debe aumentar en 0.001 en cada división.

Utilizando una base de datos “cu.summary” de precios de automoviles, paises en donde han sido fabricados, confiabilidad, consumo de gasolina, y gama del vehiculo, realizamos una estimación con un arbol de regresión.

# Regression Tree Example
library(rpart)

# grow tree
fit <- rpart(Mileage~Price + Country + Reliability + Type,method="anova", data=cu.summary)

printcp(fit) # display the results

## 
## Regression tree:
## rpart(formula = Mileage ~ Price + Country + Reliability + Type, 
##     data = cu.summary, method = "anova")
## 
## Variables actually used in tree construction:
## [1] Price Type 
## 
## Root node error: 1354.5833/60 = 22.576389
## 
## n=60 (57 observations deleted due to missingness)
## 
##            CP nsplit  rel error     xerror        xstd
## 1 0.622885266      0 1.00000000 1.03031850 0.175417788
## 2 0.132060610      1 0.37711473 0.52551930 0.100727441
## 3 0.025440940      2 0.24505412 0.36722405 0.080138910
## 4 0.011603894      3 0.21961318 0.38568677 0.090458117
## 5 0.010000000      4 0.20800929 0.40256940 0.089130973

plotcp(fit) # visualize cross-validation results

summary(fit) # detailed summary of splits

## Call:
## rpart(formula = Mileage ~ Price + Country + Reliability + Type, 
##     data = cu.summary, method = "anova")
##   n=60 (57 observations deleted due to missingness)
## 
##              CP nsplit    rel error       xerror          xstd
## 1 0.62288526607      0 1.0000000000 1.0303184965 0.17541778783
## 2 0.13206061011      1 0.3771147339 0.5255192987 0.10072744084
## 3 0.02544093971      2 0.2450541238 0.3672240504 0.08013890951
## 4 0.01160389411      3 0.2196131841 0.3856867743 0.09045811723
## 5 0.01000000000      4 0.2080092900 0.4025694027 0.08913097265
## 
## Variable importance
##   Price    Type Country 
##      48      42      10 
## 
## Node number 1: 60 observations,    complexity param=0.6228852661
##   mean=24.58333333, MSE=22.57638889 
##   left son=2 (48 obs) right son=3 (12 obs)
##   Primary splits:
##       Price       < 9446.5  to the right, improve=0.6228852661, (0 missing)
##       Type        splits as  LLLRLL,      improve=0.5044405322, (0 missing)
##       Reliability splits as  LLLRR,       improve=0.1263005365, (11 missing)
##       Country     splits as  --LRLRRRLL,  improve=0.1243525220, (0 missing)
##   Surrogate splits:
##       Type    splits as  LLLRLL,     agree=0.950, adj=0.750, (0 split)
##       Country splits as  --LLLLRRLL, agree=0.833, adj=0.167, (0 split)
## 
## Node number 2: 48 observations,    complexity param=0.1320606101
##   mean=22.70833333, MSE=8.498263889 
##   left son=4 (23 obs) right son=5 (25 obs)
##   Primary splits:
##       Type        splits as  RLLRRL,      improve=0.43853834880, (0 missing)
##       Price       < 12154.5 to the right, improve=0.25748498350, (0 missing)
##       Country     splits as  --RRLRL-LL,  improve=0.13345695140, (0 missing)
##       Reliability splits as  LLLRR,       improve=0.01637085564, (10 missing)
##   Surrogate splits:
##       Price   < 12215.5 to the right, agree=0.812, adj=0.609, (0 split)
##       Country splits as  --RRLRL-RL,  agree=0.646, adj=0.261, (0 split)
## 
## Node number 3: 12 observations
##   mean=32.08333333, MSE=8.576388889 
## 
## Node number 4: 23 observations,    complexity param=0.02544093971
##   mean=20.69565217, MSE=2.907372401 
##   left son=8 (10 obs) right son=9 (13 obs)
##   Primary splits:
##       Type    splits as  -LR--L,      improve=0.515359607900, (0 missing)
##       Price   < 14962   to the left,  improve=0.131259377800, (0 missing)
##       Country splits as  ----L-R--R,  improve=0.007022106632, (0 missing)
##   Surrogate splits:
##       Price < 13572   to the right, agree=0.609, adj=0.1, (0 split)
## 
## Node number 5: 25 observations,    complexity param=0.01160389411
##   mean=24.56, MSE=6.4864 
##   left son=10 (14 obs) right son=11 (11 obs)
##   Primary splits:
##       Price       < 11484.5 to the right, improve=0.09693168203, (0 missing)
##       Reliability splits as  LLRRR,       improve=0.07767167054, (4 missing)
##       Type        splits as  L--RR-,      improve=0.04209833909, (0 missing)
##       Country     splits as  --LRRR--LL,  improve=0.02201687475, (0 missing)
##   Surrogate splits:
##       Country splits as  --LLLL--LR, agree=0.80, adj=0.545, (0 split)
##       Type    splits as  L--RL-,     agree=0.64, adj=0.182, (0 split)
## 
## Node number 8: 10 observations
##   mean=19.3, MSE=2.21 
## 
## Node number 9: 13 observations
##   mean=21.76923077, MSE=0.7928994083 
## 
## Node number 10: 14 observations
##   mean=23.85714286, MSE=7.693877551 
## 
## Node number 11: 11 observations
##   mean=25.45454545, MSE=3.520661157

# create additional plots
par(mfrow=c(1,2)) # two plots on one page
rsq.rpart(fit) # visualize cross-validation results  

## 
## Regression tree:
## rpart(formula = Mileage ~ Price + Country + Reliability + Type, 
##     data = cu.summary, method = "anova")
## 
## Variables actually used in tree construction:
## [1] Price Type 
## 
## Root node error: 1354.5833/60 = 22.576389
## 
## n=60 (57 observations deleted due to missingness)
## 
##            CP nsplit  rel error     xerror        xstd
## 1 0.622885266      0 1.00000000 1.03031850 0.175417788
## 2 0.132060610      1 0.37711473 0.52551930 0.100727441
## 3 0.025440940      2 0.24505412 0.36722405 0.080138910
## 4 0.011603894      3 0.21961318 0.38568677 0.090458117
## 5 0.010000000      4 0.20800929 0.40256940 0.089130973

# plot tree
plot(fit, uniform=TRUE,
   main="Regression Tree for Mileage ")
text(fit, use.n=TRUE, all=TRUE, cex=.8)