ANALISIS ESPECTRAL

Series estacionaria.

Sea un conjunto de observaciones de una variable aleatoria $x$, en distintos momentos del tiempo. Cada observación de $x$, es considerada como una realización de un proceso estocástico ergódico, debido a que solo disponemos de una realización del proceso estocástico que ha generado la serie de datos, dada la imposibilidad de observar distintas realizaciones de $x(t)$ a lo largo de un periodo de tiempo.

Un proceso estocástico es estacionario en sentido estricto, cuando para todo $n>0$ la función de distribución conjunta de:

\[ F(x_{t+1},x_{t+2},...,x_{t+n})=F(x_{t+1+k},x_{t+2+k},...,x_{t+n+k}), \vee k\]

Es decir; la función de distribución conjunta es independiente de $t$ e invariante ante translaciones de tiempo.

En un sentido amplio, para que un proceso sea estacionario, es suficiente que su esperanza y su función de autocovarianza sea independiente de $t$.

Es decir,

\[E(x_t)=E(x_{t+k}), \vee k\]

Si $x_t$ es un proceso es estacionario en media entonces $\hat \mu_t=\frac 1 n \sum_{t=1}^n x_t$ es un estimador insesgado y consistente de $E(x_t)$.

Si $x_t$ es un proceso es estacionario en covarianza, se cumple la siguiente igualdad \[\gamma (t,\tau)=E([x_t-E(x_t)][x_{t+\tau}-E(x_{t+\tau})])=\gamma(\tau)\]

que significa que la función de autocovarianza no depende de t,$\gamma(\tau)=\gamma(-\tau)$, y el estimador de $\gamma(\tau)$ viene dado por

\[C(K)= \frac 1 {n-k} \sum_{t=1}^{n-k} ([x_t-\hat \mu][x_{t+k}-\hat \mu])\]

La varianza, $\gamma=0$, se estimaría a partir de:

\[C(0)= \frac 1 {n} \sum_{t=1}^{n} ([x_t-\hat \mu][x_{t}-\hat \mu])\].

Ejemplo 1

Generamos una serie aleatoria de 100 datos, con media 0 y desviación típica 1, que representamos en la figura siguiente:

x=rnorm(100,0,1)
plot(x,type='l')

Se comprueba que se trata de una serie estacionaria en media; ya que cualquier promedio que calculemos con dichos datos dará un resultado cercano a cero:

data.frame("media 1 a 25"= mean(x[1:25]),"media 26 a 50"= mean(x[26:50]),"media 51 a 75"= mean(x[51:75]),"media 76 a 100"= mean(x[76:100]))

##   media.1.a.25 media.26.a.50 media.51.a.75 media.76.a.100
## 1    0.4681097    0.01782575    -0.1115104      0.1603393

Dado que el valor esperado es cero, el estimador de la función de autocovarianza será: $C(K)= \frac 1 {n-k} \sum_{t=1}^{n-k} ([x_t][x_{t+k}])$

Calculado para diferentes valores de $k$

autocov= function(x,k) {sum(x[1:(100-k)]*x[(1+k):100])/(100-k)}
K=seq(1:10)
C_K=autocov(x,1)
for (i in 2 :10) {C=autocov(x,i)
C_K=c(C_K,C)}
data.frame(K,C_K)

##     K         C_K
## 1   1 -0.04971559
## 2   2 -0.08863767
## 3   3  0.17653272
## 4   4 -0.12609278
## 5   5  0.02457966
## 6   6  0.10952743
## 7   7 -0.04208893
## 8   8  0.15012228
## 9   9  0.06056424
## 10 10 -0.17441760

plot(K,C_K,type="h")

Como se puede apreciar el valor de C(K) es independiente de $K$, tanto en su valor como en su signo. De esta forma que el proceso aleatorio que ha generado nuestros datos es estacionario.

Si generamos a partir de estos datos una serie del tipo: $y_t=0.5+y_{t-1}$:

x0=0.5 #punto inicial 0.5
T=100
x=c(x0,rnorm(T-1))
y=cumsum(x)
plot(y,type="l")

Obtenemos una serie que no es estacionaria; ya que los promedios que obtenemos son diferentes:

data.frame("media 1 a 25"= mean(y[1:25]),"media 26 a 50"= mean(y[26:50]),"media 51 a 75"= mean(y[51:75]),"media 76 a 100"= mean(y[76:100]))

##   media.1.a.25 media.26.a.50 media.51.a.75 media.76.a.100
## 1    -2.216075     -7.077752     -8.841703       -3.99354

Y una función $C(K)= \frac 1 {n-k} \sum_{t=1}^{n-k} ([x_t-\hat \mu][x_{t+k}-\hat \mu])$ dependiente del tiempo:

autocov= function(y,k) {sum(y[1:(100-k)]*y[(1+k):100])/(100-k)}
K=seq(1:10)
C_K=autocov(y,1)
for (i in 2 :10) {C=autocov(y,i)
C_K=c(C_K,C)}
plot(K,C_K,type="h")

Análisis espectral

Una serie temporal $X_t$,es decir, una realización de un proceso estocástico puede escribirse como suma de senos y cosenos de la forma siguiente:

\[X_t=\eta+\sum_{n=1}^N\left[a_n cos(w_n t) +b_n sin(w_nt)\right] \]

donde $\eta$ es la media de la serie, $w_n=\frac {2 \pi n}{T}$ son las frecuencias naturales, cuando $T$ es par $N=\frac T 2$, y cuando $T$ es impar entonces $N=\frac {T-1} 2$.

Los coeficientes $a_n$ y $b_j$ son su amplitud y $t$ es un indice de tiempo que va de 1 a $N$.

Los coeficientes $a_n$ y $b_j$, se pueden obtener como:

$a_n=\frac 2 T \sum_{t=1}^{T} (X_t-\eta) cos (w_n t)$ y $b_n=\frac 2 T \sum_{t=1}^{T} (X_t-\eta) sin (w_n t)$

Considerando que la mayor frecuencia que se puede observar en la serie $X_t$ es la que corresponde a medio ciclo por unidad de tiempo (frecuencia de Nyquist). Una sucesiónn se puede representar como un conjunto de pares $(t, X_t)$ o como un conjunto de pares $(w_n, a^2_n + b^2_n)$. La primera representación se llama en el dominio del tiempo y la segunda en el dominio de la frecuencia. A la dinámica de las altas frecuencias (los valores mas altos de $w_n$ ) corresponden a los ciclos cortos en tanto que la dinámica de la bajas frecuencias (pequeños valores de $w_n$) van a corresponder con los ciclos largos.

El análisis espectral se utiliza 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 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.

Un armónico de frecuencia $w$ es por tanto una función de la forma:

\[a_{w}\cos(wt) +b_{w}\sin(wt)\] que da lugar a una función periódica de periodo $\frac {2\pi}{w}$.

En el análisis armónico, las series temporales no son consideradas funciones continuas como tal, sino que se obtienen a partir de una suma de $n$ ciclos con una amplitud y un periodo determinado.

En el análisis armónioco $a_n$ y $b_n$ son variables aleatorias con:

\[E(a_n)=E(b_n)=0\] \[E(a_ma_n)=E(b_mb_n)=\begin{cases} \sigma^2_n & \text{ si } m=n \\ 0 & \text{ si } n \neq m \end{cases}\] \[E(a_n,b_m)=0, \forall n m \]

En este tipo de procesos la función de autocovarianza $\gamma(\tau)$ se obtiene:

\[\gamma(\tau)=\sum_{n=1}^N \sigma^2_n\cos(w_n\tau) \] donde $\sigma^2_n$ es la varianza del armónico j-th. De manera que en $\gamma(0)=\sum_{n=1}^N \sigma^2_n$ se muestra que la varianza total del proceso es la suma de las varianzas de cada armónico.

Una manera practica de transformar una serie temporal 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:

\[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.\]

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

A continuación se presentan varias funciones R para calcular la matriz ortogonal $W$, para pasar del dominio del tiempo al dominio de la frecuencia, y viceversa:

Función MW(a)

Obtiene la matriz ortogonal, $W$, sugerida por Harvey (1978) para una serie de

MW <- function(n) {
# Author: Francisco Parra Rodr?guez 
# 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 (1:n)
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(a)

Transforma los datos del dominio del tiempo al dominio de la frecuencia pre-multiplicandolos por la matriz ortogonal,$W$, sugerida por Harvey (1978).

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

Función gdt(a)

Transforma los datos del dominio de frecuencias al dominio del tiempo pre-multiplicandolos por la matriz ortogonal, A, sugerida por Harvey (1978).

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

ESTIMACIÓN DEL PERIODOGRAMA DE UNAS SERIE ARMONICA

Serie de Fourier

Una serie de Fourier es una serie infinita que converge puntualmente a una función continua y periódica.

\[f(t)=\frac 1 2 a_0 + \sum_{n=1}^{\infty} a_{n}\cos(nw_0t) +b_{n}\sin(nw_0t)\] donde $w_0=\frac {2\pi}{T}$ se denomina frecuencia fundamental, y $a_n$ y $b_n$ se les denomina coeficientes de fourier, y la gama total de frecuencias es el Ancho de Banda de la señal.

Los coeficientes de una serie de fourier pueden calcularse gracias a la ortogonalidad de las funciones seno y coseno.

Una manera alternativa de presentar una la serie de Fourier es:

\[f(t)=c_0 + \sum_{n=1}^{\infty} c_{n}\cos(nw_0t-\theta_n)\]

siendo: $c_0=\frac {a_0}{2}$, $c_n=\sqrt {a^2_n+b^2_n}$ y $\theta_n= \arctan \frac {a_n}{b_n}$.

Partiendo de la expresión: \[a_{n}\cos(nw_0t) + b_{n}\sin(nw_0t)\] operando:

\[\sqrt {a^2_n+b^2_n}[\frac {a_{n}}{\sqrt {a^2_n+b^2_n}}\cos(nw_0t) + \frac{b_{n}}{\sqrt {a^2_n+b^2_n}}\sin(nw_0t)]\] haciendo:

\[\begin{cases} \ \frac {a_{n}}{\sqrt {a^2_n+b^2_n}}=\cos(\theta_n) \\ \frac {b_{n}}{\sqrt {a^2_n+b^2_n}}=\sin(\theta_n) \end{cases}\]

\[\theta_n= \arctan \frac {a_n}{b_n}\]

la suma se expresa solo en función del coseno:

\[ c_n[\cos(\theta_n)\cos(nw_0t) + \sin(\theta_n)\sin(nw_0t)]=c_{n}\cos(nw_0t-\theta_n)\]

Ortogonalidad

Se dicen que las funciones del conjunto $\{f_k(t)\}$ son ortogonales en el intervalo $a<t<b$, si dos fuenciones cualesquiera, $f_m(t)$ y $f_n(t)$, de dicho conjunto, cumplen:

\[\int^b_a f_m(t)f_n(t)dt=\begin{cases} & \ 0 \ {si} \ n \neq m \\ & \ r_n \ {si} \ n = m \end{cases}\]

Las funciones $sen(t)$ y $cos(t)$ son ortogonales en el intervalo $-\pi<t<\pi$, ya que :

\[\int^\pi_{-\pi} sen(t) cons (t) dt =\frac {sen^2(t)}{2}\left.[\begin{matrix} \pi\\ -\pi \end{matrix}\right.=0\]

Las funciones del conjunto:

\[\{1,cos(w_0 t),cos(2w_0t),cos(3w_0t),...,sen(w_0 t),sen(2w_0t),sen(3w_0t),...\}\] donde $w_0=\frac{2\pi}{T}$ son ortogonales en el inervalo $-\frac T 2<t<\frac T 2$, se verifica probandolo a pares:

$f_n(t)=1$ y $f_m(t)=cos(m w_0 t)$

\[\int^{\frac T 2}_{\frac {-T} 2} 1 cos (m w_0 t) dt =\frac {sen( m w_0 t)}{m w_0}\left.[\begin{matrix} \frac T 2 \\ \frac {-T} 2 \end{matrix}\right.= \frac {2sen(mw_0\frac T 2)}{mw_0}- \frac {2sen(mw_0\frac T 2)}{mw_0}=0\]

$f_n(t)=1$ y $f_m(t)=sen(m w_0 t)$

\[\int^{\frac T 2}_{\frac {-T} 2} 1 sen (m w_0 t) dt =\frac {-cos( m w_0 t)}{m w_0}\left.[\begin{matrix} \frac T 2 \\ \frac {-T} 2 \end{matrix}\right.= \frac {-1}{mw_0}[cos(mw_0 \frac T 2)-cos(mw_0 \frac T 2)]=0\]

$f_n(t)=cos(n w_0 t)$ y $f_m(t)=cos(m w_0 t)$

\[\int^{\frac T 2}_{\frac {-T} 2} cos(n w_0 t) cos (m w_0 t) dt =\begin{cases} & \ 0 \ {si} \ n \neq m \\ & \ \frac T 2 \ {si} \ n = m \neq 0 \end{cases}\]

$f_n(t)=sen(n w_0 t)$ y $f_m(t)=sen(m w_0 t)$

\[\int^{\frac T 2}_{\frac {-T} 2} sen(n w_0 t) sen (m w_0 t) dt =\begin{cases} & \ 0 \ {si} \ n \neq m \\ & \ \frac T 2 \ {si} \ n = m \neq 0 \end{cases}\] d) $f_n(t)=sen(n w_0 t)$ y $f_m(t)=sen(m w_0 t)$

\[\int^{\frac T 2}_{\frac {-T} 2} sen(n w_0 t) cos (m w_0 t) dt = 0\] para cualquier $m$ y $n$.

Cálculo de los coeficientes Fourier

Los coeficientes de fourier se calculan multiplicando $f(t)$ por $cos(mw_0 t)$ e integrando de –T/2 a T/2:

\[\int^{\frac T 2}_{\frac {-T} 2} f(t) cos (m w_0 t) dt = \frac 1 2 a_0 \int^{\frac T 2}_{\frac {-T} 2} cos (m w_0 t) dt + \sum_{n=1}^{\infty} a_n \int^{\frac T 2}_{\frac {-T} 2} cos (n w_0 t)cos (m w_0 t) dt+\sum_{n=1}^{\infty} b_n \int^{\frac T 2}_{\frac {-T} 2} sen (n w_0 t)cos (m w_0 t) dt\] que dada la ortogonalidad de las funciones de seno y coseno implica que:

\[a_0=\frac 2 T \int^{\frac T 2}_{\frac {-T} 2} f(t) dt\]

\[a_m=\frac 2 T \int^{\frac T 2}_{\frac {-T} 2} f(t) cos (m w_0 t) dt\], $m=1,2,3...$

\[b_m=\frac 2 T \int^{\frac T 2}_{\frac {-T} 2} f(t) sen (m w_0 t) dt\]

Forma compleja de la serie de Fourier

Consideremos la serie de Fourier para una función periódica $f(t)$, con periodo $T=\frac{2\pi}{w_0}$:

\[f(t)=\frac 1 2 a_0 + \sum_{n=1}^{\infty} a_{n}\cos(nw_0t) +b_{n}\sin(nw_0t)\] Es posible obtener una forma alternativa usando las fórmulas de Euler:

\[cos(nw_0t)=\frac 1 2 (e^{inw_0t}+e^{-inw_0t})\]

\[sin(nw_0t)=\frac 1 {2i} (e^{inw_0t}-e^{-inw_0t})\] sustituyendo:

\[f(t)=\frac 1 2 a_0 + \sum_{n=1}^{\infty} a_{n} \frac 1 2 (e^{inw_0t}+e^{-inw_0t}) + b_{n}\frac 1 {2i} (e^{inw_0t}-e^{-inw_0t})\] dado que $\frac 1 {i}=-i$:

\[f(t)=\frac 1 2 a_0 + \sum_{n=1}^{\infty} [\frac 1 2 (a_n-ib_n)e^{inw_0t}+\frac 1 {2} (a_n+ib_n)e^{-inw_0t}]\] definiendo:

$c_0=\frac 1 2 a_0$, $c_n=\frac 1 2 (a_n-ib_n)$ y $c_{-n}=\frac 1 2 (a_n+ib_n)$

quedaría como:

\[f(t)= \sum_{-\infty}^{\infty} c_ne^{inw_0t}\]

expresión que se conoce como forma compleja de fourier.

Y sus coeficientes $c_n$ pueden obtenerse a partir de los coeficientes $a_n$ y $b_n$, o bien a partir de :

\[c_n= \frac 1 T \int_{0}^{T}f(t)e^{inw_0t} dt\] ## Transformada de Fourier.

La Transformada de Fourier, $F(w)$ , se define para una función continua de variable real, $f(t)$ , mediante la siguiente formula:

\[F(w)= \int_{-\infty}^{\infty}f(t)e^{2\pi iwt} dt\] siendo $i=\sqrt {-1}$, $e^{2\pi iwt}=cos(2\pi wt)+isin(2\pi wt)$ y $w$ una variable que representa las distintas frecuencias.

La Transformada de Fourier es una función compleja con una parte real y otra parte imaginaria, es decir:

\[F(w)=R(w)+I(w)\]

donde $R(w)$ es la parte real e $I(w)$ es la parte imaginaria.

La representación gráfica de la función de magnitud $F(w)$ se le denomina Espectro de Fourier y se expresa en términos del modulo del número complejo:

\[F(w)=\sqrt{R^2(w)+I^2(w)}\]

y al cuadrado de dicha función $F(w)^2$ se le denomina Espectro de potencias o Espectro de Amplitud.

El gráfico de los módulos al cuadrado frente a la frecuencia es el periodograma o espectro empírico de la sucesión $f(x)$.

El periodograma recoge la contribución que tiene cada armónico a la hora de explicar la varianza de cada serie; y cada armónico esta caracterizado por la frecuencia en que tienen lugar los ciclos. Los ciclos que tienen un elevado periodo (desde que tiene lugar un máximo al siguiente máximo) tendrán una baja frecuencia y viceversa.

Se denomina esprectro de fase a:

\[\phi(w)=artg \left [\frac {-I(w)}{R(w)} \right ]\]

Cálculo del periodograma.

Consideremos la serie temporal $X_t$ de la que disponemos de un conjunto discreto y finito de observaciones T observaciones, generadas por un proceso aleatorio $x(t)$ como el descrito en el tema 1. Dado que se busca una representación de $X_t$ que correspondiente a $T$ observaciones, ajustamos los datos a un polígono trigonométrico que se asemeje a una serie de fourier, escogiendo $w_n$ como: $w_n=\frac{2\pi n}{T}$.

Es decir,

\[X_t=\frac 1 2 a_0 + \sum_{n=1}^{N} a_{n}\cos(\frac{2\pi n t}{T}) +b_{n}\sin(\frac{2\pi n t}{T})\]

\[x_t= (X_t -\frac 1 2 a_0) = (X_t - \hat \mu) = \sum_{n=1}^{N} a_{n}\cos(\frac{2\pi n t}{T}) +b_{n}\sin(\frac{2\pi n t}{T})\]

Siendo $a_0=-\frac 2 T \sum_{t=1}^{T} X_t$

La forma habitual de obtener el periodograma, es estimar por mínimos cuadrados ordinarios los coeficientes $a_n$ y $b_n$ para cada $N=\frac T 2$ armónico, si el número de observaciones $T$ es par, para cada $N=\frac {T-1} 2$ si el numero de observaciones es impar, en un modelo especificado de la siguiente forma:

$x_t = a cos(w) t + b sin (w) t + v_t$

en la que $x_t$ sería la serie armónica, $w=w_n=\frac{2 \pi n}{T}$, $T$ es el tamaño de la serie y coincide con el periodo de mayor ciclo que es posible estimar con el tamaño de la serie, $p$ indica el orden del armónico de los $\frac{T}{2}$ ciclos, $v_t$ es un residuo no explicado al que se puede considerar irrelevante (caso determinístico) o que verifica las propiedades clásicas de la perturbación de los modelos econométricos.

El periodograma o estimador del espectro se obtendría entonces a partir de la representación de $I(w_i)=\frac {T(a^2_n+b^2_n)}{4\pi}$ frente a los $n$ armónicos, en tanto que la contribución a la varianza por cada armónico sería: $\frac {(a^2_n+b^2_n)}{2}$.

Si una serie temporal de ciclos empíricos presenta en su periodograma unos pocos ciclos que explican un porcentaje significativo de su varianza, se puede obtener el ciclo teórico de dicha serie temporal a partir de los armónicos correspondientes a dichos ciclos.

Teorema de Paserval

Sea $f(x)$ una función continua en el intervalo $[-\pi ,\pi ]$ de periodo $2\pi$, con desarrollo de Fourier:

\[f(x)=\sum_{x=-\infty}^{\infty} c_n e^{inx}\] donde, $c_n$ han sido obtenidos a partir de los coeficientes $a_n$ y $b_n$.

Entonces se verifica que:

\[\sum_{n=-\infty}^{\infty} \mid c_n \mid ^2 = \frac {1}{2\pi}\int_{-\pi}^{\pi} \mid f(x) \mid^2 dx\]

Particularizando a la serie función periódica $f(t)$ , con periodo $T=\frac{2\pi}{w_o}$:

\[f(t)=\frac 1 2 a_0 + \sum_{n=1}^{\infty} a_{n}\cos(w_0nt) +b_{n}\sin(w_0nt)\]

La identidad de Paaserval quedaría:

\[ \frac {1}{2\pi}\int_{-\pi}^{\pi} [f(t)]^2=\frac 1 2 a_0^2 + \sum_{n=-\infty}^{\infty} (a_n^2+b_n^2)\] Las series temporales no son consideradas funciones continuas como tal, sino muestras de señales continuas tomadas a una misma distancia temporal a partir de un valor inicial $X_0$ y siendo $T$ el tamaño de la serie. De acuerdo a lo anterior; en la función periódica $f (t)$ la potencia promedio está dada por:

\[ \frac {1}{T}\int_{-\frac T 2}^{\frac T 2} [f(t)]^2=\frac 1 4 a_0^2 + \frac 1 2 \sum_{n=1}^{\frac T 2} (a_n^2+b_n^2)\] que muestra así que el periodograma estudia de hecho la distribución de la varianza o potencia de la serie en función de los diversos armónicos:

\[ \sigma^2= a_{\frac T 2}^2+ \frac 1 2 \sum_{n=t}^{\frac T 2} (a_n^2+b_n^2) \]

Test sobre el periodograma

Una forma de contrastar la existencia de algún ciclo en el periodograma de una serie temporal es el test de Fisher, estadístico $g$ (Fisher, 1929) o relación entre la mayor varianza asociada a una determinada frecuencia ($w_i$) y la varianza total de la serie:

\[ g=\frac{max\mid w_n \mid^2}{\frac T 2 \sum_{n=1}^{\frac T 2} \mid w_n \mid^2 }\]

Para probar la significación del periodo p se contrasta el estadístico $g$ contra la $z$ de una distribución normal (0,1), siendo la regla de decisión rechazar la hipótesis nula sobre un componente periódico en $X_t$ si la $g$ calculada excede de la $z$ en un nivel de significación del $100\alpha%$.

La manera habitual de contrastar la existencia de algún ciclo en el periodograma de una serie temporal a través del estadístico es calculando:

\[G=\frac {maxS^2}{2S^2}\]

El ciclo es significativo si el valor $G$ de esta relación es igual al valor crítico calculado según la siguiente fórmula:

\[G_c=1-e^{\frac{ln(p)-ln(T)}{T-1}}\] Siendo $ln(p)$ el logaritmo neperiano del nivel de probabilidad elegido y $T$ el número total de datos de la serie (en series de más de 30 datos).

Una prueba para estudiar la dependencia serial (Durbin; 1969) en series de observaciones estacionarias ($x_0 ,..., x_T$) se realiza sobre la grafica del periodograma acumulado:

\[s_j=\frac {\sum_{n=1}^{j}p_n}{\sum_{n=1}^{N}p_n}\]

donde $n=1,...,N$ y $p_n$ el periodograma ordinario:

\[p_n=\frac 2 T \left | {\sum_{t=1}^{T}x_t e ^{\frac { i2\pi n t}{T}}} \right | ^2\]

Se presupone que cuando $X_t$ esta independientemente y normalmente distribuida, $s_j$ se distribuye igual que el orden estadístico de $N-1$ muestras independientes de la distribución uniforme (0,1). Bartlett’s (1954,1966 p 361) sugiere para probar la independencia serial; probar la máxima discrepancia entre $s_j$ y su expectativa, ie. $\frac j N$. Para una probar un exceso de bajas frecuencias relativas frente a altas frecuencias que equivaldría a la expectativa de presencia de correlación serial positiva este enfoque conduce al estadístico:

\[ C^{+}=max(s_j - \frac j N)\]

Por el contrario un test contra excesos de variaciones de alta frecuencia el estadístico apropiado es:

\[ C^{-}=max(\frac j N- s_j)\]

El estadístico correspondiente a las dos partes de la prueba sería:

\[ C=max(C^+,C^-)\] Los valores críticos para estos estadísticos están dado en la Tabla siguiente:

Ejemplo 2

Partimos de una serie temporal generada a partir de un paseo aleatorio o random walk:$X_t = 0.5+ X_{t-1} + u_t$.

La serie $X_t$ presenta una tendencia estocástica, y vamos a descomponerla utilizando un modelo armónico, partiendo de una representación de la tendencia ó movimiento relevante de la serie temporal obtenida a partir de una tendencia cuadrática, $k=\frac T 2$ ciclos armónicos y un residuo aleatorio $v_t$:

\[X_t=a+bt+ct^2 + \sum_{n=1}^{N} a_{n}\cos(w_0nt) +b_{n}\sin(w_0nt)+v_t\]

de manera que:

\[x_t= X_t -a+bt+ct^2= \sum_{n=1}^{N} a_{n}\cos(w_0nt) +b_{p}\sin(w_0nt)+v_t\]

set.seed(50)
n=100
w = rnorm(n) 
x = 0.5+cumsum (w)
plot(x,t="l")

t=seq(1:n)
tendencia=lm(x~t+t^2)$fitted
plot(t,x,type="l")
lines(t,tendencia,col="red")

El armónico de periodo $1$ (el que tinen un mínimo y un máximo en el espacio temporal de la serie) se elabora a partir de $\frac {cos(2 \pi t)}{100}$ y $\frac {sin(2 \pi t)}{100}$ para $t=1,...,100$. La representación gráfica de ambas series aparece en la figura siguiente:

armo1.1=cos(2*pi*t/n)
armo1.2=sin(2*pi*t/n)

plot(t,armo1.1,type="l")
lines(t,armo1.2,col="red")

La regresión minimo cuadrática entre ambas series y la serie libre de tendencia ($x_t$); ofrece el siguiente resultado que se acompaña con su representación gráfica:

lm((x-tendencia)~0+armo1.1+armo1.2)

## 
## Call:
## lm(formula = (x - tendencia) ~ 0 + armo1.1 + armo1.2)
## 
## Coefficients:
## armo1.1  armo1.2  
## -0.8196   1.6038

armo1=lm((x-tendencia)~0+armo1.1+armo1.2)$fitted
plot(t,armo1,type="l")

Este proceso repetido para los 50 periodos permite obtener los coeficientes con los que elaborar el peridograma y obtener la contribución de cada armónico a la varianza de la serie:

armo=function(a,b) {
t=seq(1:length(a))
p=length(a)/b
tendencia=lm(a~t+t^2)$fitted
armo.1=cos(2*pi*t/p)
armo.2=sin(2*pi*t/p)
reg=lm((a-tendencia)~0+cos(2*pi*t/p)+sin(2*pi*t/p))
data.frame(frecuencia=b,periodo=p,a=coef(reg)[1],b=coef(reg)[2])
}
res=armo(x,1)
for(i in 2:20) {
  res=rbind(res,armo(x,i))} 

periodogr=data.frame(frecuencia=res$frecuencia,periodo=res$period,a=res$a,b=res$b,Pr=n*(res$a^2+res$b^2)/(4*pi),contribución=(res$a^2+res$b^2)/(2))

periodogr

##    frecuencia    periodo           a            b          Pr contribución
## 1           1 100.000000 -0.81955539  1.603756807 25.81259964 1.6218534683
## 2           2  50.000000  2.32667743 -1.554643932 62.31191040 3.9151727990
## 3           3  33.333333  0.84435020 -0.111193759  5.77168491 0.3626456583
## 4           4  25.000000  0.07134080 -0.699086441  3.92962596 0.2469056812
## 5           5  20.000000  0.04605747 -0.325094722  0.85790777 0.0539039347
## 6           6  16.666667  0.35212968 -0.872266756  7.04136964 0.4424223027
## 7           7  14.285714  0.72134651 -0.459005208  5.81732456 0.3655132820
## 8           8  12.500000  0.03514291  0.017685607  0.01231704 0.0007739023
## 9           9  11.111111 -0.31723145 -0.523382452  2.98069342 0.1872824908
## 10         10  10.000000  0.08707847 -0.420452596  1.46711450 0.0921815225
## 11         11   9.090909 -0.13988519  0.007170838  0.15612532 0.0098096434
## 12         12   8.333333 -0.03667392 -0.073060020  0.05317958 0.0033413716
## 13         13   7.692308  0.04226489  0.204201892  0.34604051 0.0217423667
## 14         14   7.142857  0.13214521 -0.240579621  0.59954392 0.0376704558
## 15         15   6.666667  0.16110274 -0.083340076  0.26180717 0.0164498298
## 16         16   6.250000  0.13108038  0.016066559  0.13878471 0.0087201003
## 17         17   5.882353 -0.02683302 -0.267789409  0.57638900 0.0362155892
## 18         18   5.555556 -0.06243762 -0.038746417  0.04296977 0.0026998705
## 19         19   5.263158  0.18952061 -0.264605838  0.84299846 0.0529671552
## 20         20   5.000000 -0.14339838  0.131332526  0.30089298 0.0189056633

plot(periodogr$Pr,type="l",col="red")

Para comprobar la significación estadística del ciclo de o periodo $50$, calculamos el estadístico $G=\frac {maxS^2}{2S^2}$:

G=max(periodogr$Pr)^2/(2*periodogr$Pr[2]^2)
G

## [1] 0.5

El ciclo es significativo para un nivel de probabilidad del 95% ya que el valor del estadístico es superior al valor crítico $G_c=1-e^{\frac{ln(0,05)-ln(50)}{99}}$

Gc=1-exp((log(0.05)-log(50))/49)
Gc

## [1] 0.1314886

Calculo del periodograma a través de la Transformada Discreta de Fourier

Tomando $N$ muestras de una señal periodica $x_k = f(t_k)$ de periodo $T$ en instantes separados por intervalos regulares:

\[t_0=0, t_1=\frac T N, t_2=\frac {2T} N,...,t_k=\frac {kT} N,..., t_{N-1}=\frac {(N-1)T} N\] Cabe aproximarla mediante una combinación $g(t)$ de funciones T-periódicas conocidas que tome en dichos puntos el mismo valor que $f$. Este procedimiento se conoce como interpolación trigonométrica. Las funciones T-periódicas que se utilizan son los armónicos complejos $e^{inwt}$ con $w=\frac {2\pi} T$ y puesto que hay $N$ puntos, si el problema tiene solución única hay que combinar un total de $N$ armónicos.

La función $g(t)$ utilizada en la aproximación, toma entonces la forma general:

\[g(t)=\frac 1 N (\beta_0+\beta_1e^{iwt}+\beta_2e^{i2wt}+....+\beta_ke^{ikwt}+...+\beta_{N-1}e^{i(N-1)wt})=\frac 1 N \sum_{n=0}^{N-1}\beta_i e^{inwt}\] tal que $X_k = g(t_k)$ para cada $k=0,1,...,N-1$.

Entonces,

\[x_k=\frac 1 N \sum_{\eta=0}^{N-1}\beta_{\eta} e^{i\eta wt_k}=\frac 1 N \sum_{\eta=0}^{N-1}\beta_{\eta} e^{\frac{i\eta w2\pi}{N}}=\frac 1 N \sum_{\eta=0}^{N-1}\beta_{\eta} w^{nk}_N\]

siendo $w_N=e^{\frac{i2\pi}{N}}$ la raiz primitiva N-esima de la unidad.

En forma matricial se expresa:

\[ \left(\begin{array}{c} x_0 \\ x_1 \\ x_2 \\ . \\ x_{\eta} \\ . \\ x_{N-1} \\ \end{array} \right)= \frac{1}2\left( \begin{array}{ccccccc} 1& 1& 1& ... & 1 & ... & 1 \\ 1& w& w^2 & ... & w^{\eta} & ... & w^{N-1} \\ 1& w^2& w^4& ... & w^{2\eta} & ... & w^{2(N-1)} \\ .& .&. & ...&. & ...& . \\ 1& w^k& w^{2k}& ... & w^{\eta k} & ... & w^{(N-1)k} \\ .& .&. & ...&. & ...& . \\ \\ 1& w^{N-1}& w^{2(N-1)}& ... & w^{\eta(N-1)} & ... & w^{(N-1)(N-1)} \\ \end{array} \right) \left(\begin{array}{c} \beta_0 \\ \beta_1 \\ \beta_2 \\ . \\ \beta_{\eta} \\ . \\ \beta_{N-1} \\ \end{array} \right) \] donde $F_N=[w^{nk}]^{N-1}_{n,k=0}$ es la matriz de Fourier de orden N.

Al vector $\beta$ se le denomina transformada discreta de Fourier del vector $x$, denotandose como : $\beta = DFT(x)$.

Una forma de obtener la DFT es a través del algoritmo FFT (Fast Fourier Transform), desarrollado por diseñado por J.W. Cooley y John Tukey en 1965.

Si la función que interpolamos es una función real de periodo $T$, $x_k=g(t_k)$, donde $k = 0,1,..., N -1$, que utiliza la forma general:

\[g(t)=\sum_n a_ncos(nwt)+b_nsin(nwt) \] con $w=\frac{2\pi}{T}$ suponiendo que $M=2N$,si $\beta = DFT(x)$, entonces:

\[a_0=\frac {\beta_0}{N}; a_n=\frac {2RE(\beta_n)}{N};b_n=\frac {-2IM(\beta_n)}{N}, N=1,...,M-1;a_M=\frac {\beta_M}{N}\] y el poligono trigonométrico quedaría:

\[g(t)=a_0+\sum_{n=1}^{M-1} a_ncos(nwt)+b_nsin(nwt) + a_Mcos(Mwt) \]

Si se desea hacer una gráfica acumulativa del periodograma (B.D. Ripley), existe la función “cpgram”, en esta prueba se utilizan como bandas de frecuencia máxima y mínima las siguientes:

\[max= c+(\frac t n) ; min= -c+(\frac t n)\]

donde $c=\frac{1.358}{\sqrt {n-1}+0.12+ \frac {0.11}{\sqrt{n-1}}}$ y $n$ el número de frecuencias que se utilizan en el periodograma.

Si la representación grafica acumulativa del periodograma situada entre situada las bandas de frecuencia máxima y mínima, indicaría la no presencia de ciclos dominantes en la serie temporal.

Ejemplo 3

Utilizamos la serie ajustada a tendencia del Ejemplo 1, calculamos la tranformada de fourier con la función generica de R “fft”, la serie se obtendría a través de la inversa:

x.1=x-tendencia
z <- fft(x.1)

Para obtener el periodograma:

CF = abs(z/sqrt(100))^2
P = (4/100)*CF[1:51] # Solo se necesitan los (n/2)+1 valores de la FFT .
f=(0:50)/16 # Para crear las frecuencias armónicas de 1/100 en pasos de 0 a 0.5  
plot(f,P,type="l") # gráfica del periodograma; tipo = “l” gráficos de línea.

Se puede calcular directamente el periodograma utilizando la transformada de Fourier con las funcion de R: “spec.pgram”:

spec.pgram(x.1)

Realizamos la representación gráfica del periodograma a través de "cpgram$.

cpgram(x.1)

Cálculo del peridograma transformando la serie al dominio de la frecuencia.

Sea $x$ un vector T x 1, el modelo transformado en el dominio de la frecuencia por la matriz ortogonal $W$, esta dado por:

$\dot x= Wx$

Transformado las series del dominio del tiempo al dominio de la frecuencia, puede obtenerse un estimador consistente del periodograma de $x$.

Denominando $p_n$ el ordinal del periodograma de $x$ en la frecuencia $\lambda_n=\frac{2\pi j} T$, y $\dot x_n$ el enesimo elemento de $\dot x$ :

\[ \left\lbrace \begin{array}{ll} p_{2n}=\dot x_{2n}^{2}+\dot x_{2n+1}^{2} & \forall n = 1,...\frac{T-1}{2}\\ p_{2n}=\dot x_{2n}^{2}& \forall n = \frac{T-1}{2} \end{array} \right .\]

\[p_1=\dot x_{1}^{2}\]

Las siguientes funciones R calculan y representan gráficamente el periodograma:

Función periodograma (a)

Calcula y presenta el espectro de la serie “a”

periodograma <- function(y) {
# Author: Francisco Parra Rodríguez
# Some ideas from Gretl 
# http://econometria.wordpress.com/2013/08/21/estimation-of-time-varying-regression-coefficients/ 
cf <- gdf(y)
n <- length(y)
if (n%%2==0) {
m1 <- c(0)
m2 <- c()
for(i in 1:n){
if(i%%2==0) m1 <-c(m1,cf[i]) else m2 <-c(m2,cf[i])}
m2 <-c(m2,0) 
frecuencia <- seq(0:(n/2)) 
frecuencia <- frecuencia-1
omega <- pi*frecuencia/(n/2)
periodos <- n/frecuencia
densidad <- (m1^2+m2^2)/(4*pi)
tabla <- data.frame(omega,frecuencia, periodos,densidad)
tabla$densidad[(n/2+1)] <- 4*tabla$densidad[(n/2+1)]
data.frame(tabla[2:(n/2+1),])}
else {m1 <- c(0)
m2 <- c()
for(i in 1:(n-1)){
if(i%%2==0) m1 <-c(m1,cf[i]) else m2 <-c(m2,cf[i])}
m2 <-c(m2,cf[n]) 
frecuencia <- seq(0:((n-1)/2)) 
frecuencia <- frecuencia-1
omega <- pi*frecuencia/(n/2)
periodos <- n/frecuencia
densidad <- (m1^2+m2^2)/(4*pi)
tabla <- data.frame(omega,frecuencia, periodos,densidad)
data.frame(tabla[2:((n+1)/2),])}
}

Función gperiodograma (a)

Presenta gráficamente el espectro de la variable a

gperiodograma <- function(y) {
# Author: Francisco Parra Rodr?guez
# Some ideas from Gretl 
# http://econometria.wordpress.com/2013/08/21/estimation-of-time-varying-regression-coefficients/ 
tabla <- periodograma(y)
plot(tabla$frecuencia,tabla$densidad,
main = "Espectro", 
ylab = "densidad",
xlab="frecuencia",type = "l",
col="#ff0000")}

Las funciones R para realizar el test de Durbin sobre el periodograma son:

Función td (y,significance)

El test de Durbin esta basado en el siguiente estadistico: $s_j=\frac{\sum_{n=1} ^j p_n}{\sum_{n=1}^N p_n}$

donde $N=\frac{T}{2}$ para $T$ par y $\frac{T-1}{2}$ para $T$ 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$

En la función se utilizan las bandas de frecuencia de la función cpgram de B. D. Ripley.

td <- function(y) {
# Author: Francisco Parra Rodríguez
# Some ideas from:
#Harvey, A.C. (1978), Linear Regression in the Frequency Domain, International Economic Review, 19, 507-512.
# DURBIN, J., "Tests for Serial Correlation in Regression Analysis based on the Periodogram ofLeast-Squares Residuals," Biometrika, 56, (No. 1, 1969), 1-15.
# Venables and Ripley, "Modern Applied Statistics with S" (4th edition, 2002).
# http://econometria.wordpress.com/2013/08/21/estimation-of-time-varying-regression-coefficients/ 
per <- periodograma(y)
p <- as.numeric(per$densidad)
n <- length(p)
s <- p[1]
t <- 1:n
for(i in 2:n) {s1 <-p[i]+s[(i-1)]
s <- c(s,s1)
s2 <- s/s[n]
}
xm <- frequency(y)/2
c <- 1.358/(sqrt(n-1)+0.12+0.11/sqrt(n-1))
min <- -c+(t/length(p))
max <- c+(t/length(p))
data.frame(s2,min,max)
}

Función gtd (a)

Presenta gráficamente los resultados de la prueba de Durbin (Durbin; 1969) :

gtd <- function (y) {
S <- td(y)
plot(ts(S), plot.type="single", lty=1:3,main = "Test Durbin", 
ylab = "densidad acumulada",
xlab="frecuencia")
}

Ejemplo 4

Obtenemos el periodograma, su representación gráfica y el test de Durbin sobre el periodograma utilizando de la serie creada en del ejemplo2.

periodograma(x.1)

##         omega frecuencia   periodos     densidad
## 2  0.06283185          1 100.000000 1.290630e+01
## 3  0.12566371          2  50.000000 3.115596e+01
## 4  0.18849556          3  33.333333 2.885842e+00
## 5  0.25132741          4  25.000000 1.964813e+00
## 6  0.31415927          5  20.000000 4.289539e-01
## 7  0.37699112          6  16.666667 3.520685e+00
## 8  0.43982297          7  14.285714 2.908662e+00
## 9  0.50265482          8  12.500000 6.158519e-03
## 10 0.56548668          9  11.111111 1.490347e+00
## 11 0.62831853         10  10.000000 7.335572e-01
## 12 0.69115038         11   9.090909 7.806266e-02
## 13 0.75398224         12   8.333333 2.658979e-02
## 14 0.81681409         13   7.692308 1.730203e-01
## 15 0.87964594         14   7.142857 2.997720e-01
## 16 0.94247780         15   6.666667 1.309036e-01
## 17 1.00530965         16   6.250000 6.939235e-02
## 18 1.06814150         17   5.882353 2.881945e-01
## 19 1.13097336         18   5.555556 2.148489e-02
## 20 1.19380521         19   5.263158 4.214992e-01
## 21 1.25663706         20   5.000000 1.504465e-01
## 22 1.31946891         21   4.761905 1.798089e-01
## 23 1.38230077         22   4.545455 1.673123e-02
## 24 1.44513262         23   4.347826 5.822690e-01
## 25 1.50796447         24   4.166667 1.364430e-01
## 26 1.57079633         25   4.000000 1.932048e-01
## 27 1.63362818         26   3.846154 5.732373e-03
## 28 1.69646003         27   3.703704 8.312120e-02
## 29 1.75929189         28   3.571429 1.134812e-01
## 30 1.82212374         29   3.448276 1.975984e-02
## 31 1.88495559         30   3.333333 2.045660e-02
## 32 1.94778745         31   3.225806 1.900729e-02
## 33 2.01061930         32   3.125000 7.589156e-02
## 34 2.07345115         33   3.030303 7.107889e-02
## 35 2.13628300         34   2.941176 8.303240e-02
## 36 2.19911486         35   2.857143 1.216226e-01
## 37 2.26194671         36   2.777778 4.915373e-02
## 38 2.32477856         37   2.702703 4.859259e-02
## 39 2.38761042         38   2.631579 2.036695e-02
## 40 2.45044227         39   2.564103 5.476973e-03
## 41 2.51327412         40   2.500000 1.581200e-02
## 42 2.57610598         41   2.439024 5.910273e-02
## 43 2.63893783         42   2.380952 1.466372e-01
## 44 2.70176968         43   2.325581 7.117237e-04
## 45 2.76460154         44   2.272727 2.598291e-02
## 46 2.82743339         45   2.222222 1.473444e-01
## 47 2.89026524         46   2.173913 7.863477e-02
## 48 2.95309709         47   2.127660 7.153018e-02
## 49 3.01592895         48   2.083333 2.445113e-02
## 50 3.07876080         49   2.040816 1.296646e-01
## 51 3.14159265         50   2.000000 2.956689e-01

gperiodograma(x.1)

td(x.1)

##           s2         min       max
## 1  0.2064961 -0.17031031 0.2103103
## 2  0.7049802 -0.15031031 0.2303103
## 3  0.7511526 -0.13031031 0.2503103
## 4  0.7825889 -0.11031031 0.2703103
## 5  0.7894520 -0.09031031 0.2903103
## 6  0.8457817 -0.07031031 0.3103103
## 7  0.8923192 -0.05031031 0.3303103
## 8  0.8924178 -0.03031031 0.3503103
## 9  0.9162628 -0.01031031 0.3703103
## 10 0.9279994  0.00968969 0.3903103
## 11 0.9292484  0.02968969 0.4103103
## 12 0.9296738  0.04968969 0.4303103
## 13 0.9324421  0.06968969 0.4503103
## 14 0.9372383  0.08968969 0.4703103
## 15 0.9393327  0.10968969 0.4903103
## 16 0.9404430  0.12968969 0.5103103
## 17 0.9450540  0.14968969 0.5303103
## 18 0.9453978  0.16968969 0.5503103
## 19 0.9521416  0.18968969 0.5703103
## 20 0.9545487  0.20968969 0.5903103
## 21 0.9574256  0.22968969 0.6103103
## 22 0.9576933  0.24968969 0.6303103
## 23 0.9670093  0.26968969 0.6503103
## 24 0.9691924  0.28968969 0.6703103
## 25 0.9722836  0.30968969 0.6903103
## 26 0.9723753  0.32968969 0.7103103
## 27 0.9737052  0.34968969 0.7303103
## 28 0.9755209  0.36968969 0.7503103
## 29 0.9758370  0.38968969 0.7703103
## 30 0.9761643  0.40968969 0.7903103
## 31 0.9764684  0.42968969 0.8103103
## 32 0.9776827  0.44968969 0.8303103
## 33 0.9788199  0.46968969 0.8503103
## 34 0.9801484  0.48968969 0.8703103
## 35 0.9820943  0.50968969 0.8903103
## 36 0.9828808  0.52968969 0.9103103
## 37 0.9836582  0.54968969 0.9303103
## 38 0.9839841  0.56968969 0.9503103
## 39 0.9840717  0.58968969 0.9703103
## 40 0.9843247  0.60968969 0.9903103
## 41 0.9852703  0.62968969 1.0103103
## 42 0.9876165  0.64968969 1.0303103
## 43 0.9876278  0.66968969 1.0503103
## 44 0.9880436  0.68968969 1.0703103
## 45 0.9904010  0.70968969 1.0903103
## 46 0.9916592  0.72968969 1.1103103
## 47 0.9928036  0.74968969 1.1303103
## 48 0.9931948  0.76968969 1.1503103
## 49 0.9952694  0.78968969 1.1703103
## 50 1.0000000  0.80968969 1.1903103

gtd(x.1)

Ventanas

Hasta ahora hemos supuesto que las frecuencias eran frecuencias de Fourier y por tanto $w=w_n=\frac {2\pi n} T$, donde $n$ indica el orden del armónico. Se y se interpreta como el número de veces que un sinusoide (un armónico) de frecuencia $w_n$ ejecuta un ciclo completo en la muestra considerada. Es decir, si $n = 4$, la frecuencia asociada al armónico ($w_4=\frac {2 \pi 4} T$), determinara que este ejecute 4 ciclos completos a lo largo de $T$. A este tipo de frecuencias se denominan frecuencias de Fourier.

Si suponemos que existe un armónico que se repite cuatro veces y media,dicha frecuencia no producirá ciclos enteros en la muestra y nos encontramos con una frecuencia que no es de Fourier. Estas frecuencias originan un problema que se denomina leakage o distorsión, que determina que los pesos significativos del periodograma se repartan entre frecuencias contiguas.

Una de las maneras de solucionar el leakage consiste en aplicar transformar la serie original multiplicándola por una expresión que se denominan Data Windows o taper, y obtener el periodograma a partir de la serie transformada. Así es estimador de la función de densidad espectral puede considerarse como:

\[\hat f(w) =I(w) \omega\] DoDonde $\omega$ es la función de pesos o ventana espectral y $I(W)$ es el periodograma.

Dado de que lo que se trata es de promediar algunos valores contiguos del periodograma, podría utilizarse una media móvil de amplitud $m$:

\[ \omega_t=\left\lbrace \begin{array}{ll} \frac 1 m \forall t = 0 \pm 1\pm 2\pm ...\pm \frac {m-1} 2 \\ 0 \end{array} \right .\]

Entre las ventanas más utilizadas citar:

Tukey

\[\omega_t=1-2a+2a cos \left(\frac {\pi t} T\right)\] si $a=\frac 1 4$ se denomina ventana de Tukey-Hammond

Parcen

\[ \omega_t=\left\lbrace \begin{array}{ll} 1- 6 (\frac t T)^2 + 6 \mid \frac t T \mid^3 \forall t = 1, 2,...,\frac {T} 2 \\ 2(1 - \mid \frac t T \mid)^3 \forall t = \frac {T} 2,...,T \end{array} \right .\]

Boxcar

\[ \omega_t=\left\lbrace \begin{array}{ll} \frac 1 2 \left[ 1- cos \left(\frac {\pi \frac {t+1} 2} m \right)\right] \forall t = 1, 2,...,m \\ 1 , m \\ \frac 1 2 \left[ 1- cos \left(\frac {2\pi \frac {t+1} 2} m\right)\right] \forall t = T-m+1,...,T \end{array} \right .\]

donde $m$ es arbitrario, si bien suele elegirse un valor de $m$ tal que $\frac {2m} T$ se sitúe entre $0.1$ y $0.2$.

ANALISIS ARMONICO DE PROCESOS BIVARIANTES

Proceso bivariante

Un proceso bivariante $z(t)$ es un par formado por dos procesos univariantes, $x(t)$ y $y(t)$, donde $E[x(t)] = \mu_x(t)$ y $E[y(t)] = \mu_y(t)$.

La función de autocovarianza de $x(t)$ será:

\[\gamma_x (t,\tau) =E[(x(t)- \mu_x(t))(x(t+\tau)- \mu_x(t+\tau))]\]

en tanto que la función de autocovarianza de $y(t)$ será:

\[\gamma_y (t,\tau) =E[(y(t)- \mu_y(t))(y(t+\tau)- \mu_y(t+\tau))]\] Se denomina función de cross-varianza o covarianza cruzada a:

\[\gamma_{x,y} (t,\tau) =E[(x(t)- \mu_x(t))(y(t+\tau)- \mu_y(t+\tau))]\]

Hay que señalar que $\gamma_{x,y} (t,\tau)$ no es igual a $\gamma_{y,x} (t,\tau)$, pero existe una relación entre las dos funciones ya que:

\[\gamma_{x,y} (t,\tau)=\gamma_{y,x} (t+\tau,-\tau)\]

Señalar, por último; que la covarianza entre $x(t)$ y $y(t)$ sería $\gamma_{x,y} (t,0)$.

Si se asume la estacionariedad de $x(t)$ y $y(t)$, entonces $E[x(t)] = \mu_x$ y $E[y(t)] = \mu_y$, y la función de cross-varianza no dependerá más que del retardo $\tau$.

Suponiendo que $\mu_x=\mu_y=0$, se comprueba que $\gamma_{x,y}$, no depende mas que del retardo $\tau$, es decir que $\gamma_{x,y} (t,\tau)=\gamma_{x,y} (\tau)$:

\[\gamma_{x,y} (t,\tau) =E[x(t)y(t+\tau)]=E[x(t+s)y(t+s+\tau)],\forall s,t\]

La función de correlación cruzada se define como:

\[\rho_{xy}(\tau)=\frac{\gamma_{x,y} (\tau)}{\gamma_{x} (0)\gamma_{y} (0)}\] Cuando $t = 0$ , $\gamma_{x,y} (0)$ es la covarianza habitual y

\[\rho_{xy}(0)=\frac{\gamma_{x,y} (0)}{\gamma_{x} (0)\gamma_{y} (0)}\] sería el coeficiente de correlación de Pearson entre $x(t)$ e $y(t)$.

Los estimadores de $\gamma_{x,y} (\tau)$ y $\rho_{x,y} (\tau)$ se calculan:

\[ C_{x,y}(k) =\begin{cases} & \ \frac 1 T \sum^{T-k}_{t=1} (x(t)-\mu_x)(y(t+k)-\mu_y) \ si \ k=0,1,...,T-1 \\ & \ \frac 1 T \sum^{T-k}_{t=1} (x(t+k)-\mu_x)(y(t)-\mu_y) \ si \ k=-1,-2,...,-(T-1) \end{cases}\] \[r_{xy}(k)=\frac{C_{x,y} (k)}{\sqrt{C_{x} (0)C_{y} (0)}}\] ## Análisis armónico de un proceso bivariante.

La función de autocovarianza que obtenemos en el dominio temporal, tiene también su correspondiente representación en el dominio frecuencial; esta es el cross-espectro o espectro cruzado. Así, si partimos de dos procesos estacionarios $x(t)$ y $y(t)$, con la siguiente representación espectral:

\[x(t)=\int_{0}^{\pi}cos(wt)dU_x(w)+\int_{0}^{\pi}sin(wt)dV_x(w)\] \[y(t)=\int_{0}^{\pi}cos(wt)dU_y(w)+\int_{0}^{\pi}sin(wt)dV_y(w)\] donde $U_i(w)$ e $V_i(w)$, $i=x,y$ son procesos estocásticos con dominio definido en $(0,\pi)$, con media $0$ y de incrementos incorrelacionados.

Dado que dichos procesos son conjuntamente estacionarios en covarianza, se demuestra que:

\[E[dU_x(w)dU_y(w')]=E[dV_x(w)dV_y(w')]=E[dU_x(w)dV_y(w')]=E[dV_x(w)dU_y(w')]=0,\forall w\neq w'\] \[E[dU_x(w)dU_y(w)]=E[dV_x(w)dV_y(w)]=C(w)dw\] \[E[dU_x(w)dV_y(w)]=E[dV_x(w)dU_y(w)]=q(w)dw\]

De manera que la notación de la cross-varianza quedaría como:

\[\gamma_{x,y}(\tau)=\int_{0}^{\pi}cos(wt)C(w)dw+\int_{0}^{\pi}sin(wt)q(w)dw\] Que implica que la covarianza entre $x(t)$ e $y(t)$ sea:

\[\gamma_{x,y}(0)=\int_{0}^{\pi}C(w)dw\] El cross-espectro se formula como:

\[f_{xy}(w)=\frac 1 {\pi} \sum_{\pi=-\infty}^{\infty}\gamma_{xy}(\tau) e^{-iw\tau},0 \leq w \leq \pi \]

Dado que en general el cross-espectro es complejo; se define el cross-espectro ($C$) como la parte real de cross-espectro y el espectro de cuadratura ($Q$) como la parte imaginaria, que además coinciden con $C(w)$ y $q(w)$:

\[f_{xy}(w)=C(w)-iq(W) \]

Se deduce que:

\[C(w)=\frac 1 {\pi} \sum_{\pi=-\infty}^{\infty}\gamma_{xy}(\tau)cos(w\tau)\] y

\[q(w)=\frac 1 {\pi} \sum_{\pi=-\infty}^{\infty}\gamma_{xy}(\tau)sin(w\tau)\]

La representación trigonométrica del cross-espectro será:

\[f_{xy}(w)=\alpha_{xy}(w)e^{i\phi_{xi}(w)} \]

siendo

\[\alpha_{xy}(w)=\sqrt{C^2(w)+q^2(w)}\] que se conoce como espectro de cross-amplitud.

\[\phi_{xy}=artg \left [\frac {-q(w)}{C(w)} \right ]\] que se denomina espectro de fase.

Del cross-espectro y de la función de densidad espectral individual de las dos series $x(t)$ e $y(t)$ se obtiene la función de coherencia:

\[R(w)=\frac{C^2(w)+q^2(w)}{f_x(w) f_y(w)}\]

El cross-espectro representa la aportación a la covarianza entre $x(t)$ e $y(t)$ de sus diversos componentes armónicos. Como su interpretación no es simple, se utilizan las funciones de espectro de fase y coherencia, ya que el espectro de fase revela el desfase o retardo que en el comportamiento cíclico sigue una serie respecto a la otra; y el análisis de la función de coherencia permite identificar si la correlación que se da entre las dos series se debe a que ambas siguen un comportamiento cíclico en determinados periodos, permitiendo identificar la duración o periodo de los armónicos que dominan en ambas series a la vez y que producen una alta correlación.

La construcción del cross-espectro cuando $\tau = 0$ y $\gamma_{xy}(0)$ es la covarianza habitual, da lugar a las siguientes funciones $C(w)$ y $q(w)$ :

\[C(w)=\frac{\gamma_{xy}(0)}{2\pi} \] \[q(w)=0\]

Si $E[x(t)]=\mu_x=0$ y $E[y(t)]=\mu_y=0$, la covarianza entre $x(t)=x_t$ e $y(t)=y_t$ se reducuría a $\gamma_{xy}(0)=\sum_{t=1}^{T} x_ty_t$, y la parte real del cross-espectro se reduce a :

\[C(w)=\frac 1 {2\pi}\sum_{t=1}^{T} x_ty_t\]

Teorema de Plancharel

Sean $A(x)$ y $B(x)$ dos funciones continuas de periodo $2\pi$ cuyos desarrollos de Fourier son:

\[A(x)=\sum_{x=-\infty}^{\infty}a_ne^{inx}\] \[B(x)=\sum_{x=-\infty}^{\infty}b_ne^{inx}\]

Entonces se verifica la relación de Plancharel entre los correspondientes productos escalares:

\[\sum_{n=-\infty}^{\infty}a_nb_n= \frac {1}{2\pi}\int_{-\pi}^{\pi}A(x)B(X)dx\] Si $A(x)=B(x)$ se obtiene la identidad de Parseval:

\[\sum_{n=-\infty}^{\infty}|a_n|^2= \frac {1}{2\pi}\int_{-\pi}^{\pi}|A(x)|^2dx\]

De igual manera que la identidad de Parseval estudia la distribución de la varianza de una serie desarrollada en sus armónicos, la de Plancharel estudia la covarianza entre dos series desarrolladas en sus armónicos.

Partiendo de una serie armónica $x(t)= \sum_{n=1}^{k} a_{n}\cos(nw_0t) +b_{n}\sin(nw_0t)$ y otra $y(t)= \sum_{n=1}^{k} a_{n}^*\cos(nw_0t) +b_{n}^*\sin(nw_0t)$, en donde $k=\frac {T}{2}$ si el numero de observaciones $T$ es par, o $k=\frac {T-1}{2}$ si el numero de observaciones $T$ es impar, la expresión de igualdad de plancharel sería:

\[\frac 1 2 \sum_{n=1}^{\frac T 2} a_na_n^*+b_nb^*_n= \frac {1}{T}\int_{-\frac T 2}^{\frac T 2} x(t)y(t)dx\]

El producto escalar de x(t) e y(t):

\[\sum_{t=1}^{T} x(t) y(t) = \sum_{t=1}^{T} (\sum_{n=1}^{k} [a_{n}\cos(nw_0t) +b_{n}\sin(nw_0t) ][a_{n}^*\cos(nw_0t) +b_{n}^*\sin(nw_0t)])\] da como resultado:$\frac 1 2 \sum_{n=1}^{\frac T 2} (a_na_n^*+b_nb^*_n)$, en base a la ortogonalidad de las series de seno y coseno.

Coeficiente de correlación de Pearson

Dado que la covarianza entre las series armónicas $x(t)$ e $y(t)$ se desarrolla a partir de los coeficientes de Fourier:

\[\sigma^2_{xy}=\frac 1 2 \sum_{n=1}^{\frac T 2} a_na_n^*+b_nb^*_n\]

cabe considerar a cada expresión $\frac {a_na_n^*+b_nb^*_n} 2$ como la contribución del armónico $n$ a la formación de la covarianza, de manera que la representación de $C_{xy}(w_n)=\frac {T(a_na_n^*+b_nb^*_n)} {4\pi}$ frente a los $n$ armónicos permite apreciar las frecuencias entre las que las series $x(t)$ e $y(t)$ covarían en sentido positivo o negativo. Se puede observar que un ciclo relevante en ambas series originará un valor alto en $C_{xy}(w_n)$, en tanto que un ciclo poco relevante en alguna de las dos series dará lugar a un valor bajo en $C_{xy}(w_n)$.

En tanto que el coeficiente de correlación de pearson se obtendría a partir de:

\[\rho_{xy}(0)=\frac {\sum_{n=1}^{\frac T 2} a_na_n^*+b_nb^*_n}{\sqrt{(\sum_{n=1}^{\frac T 2} a_n^2+b_n^2})(\sum_{n=1}^{\frac T 2} a_n^{*2}+b_n^{*2})}\] Utilizando la definición alternativa de las series de fourier, tenemos que $x(t)=C_0+\sum_{n=1}^{k}C_n(cos(nw_0t-\theta_n)$ e $y(t)=C_0^*+\sum_{n=1}^{k}C_n^*(cos(nw_0t-\theta_n^*)$, en donde $C_0=\frac {a_0}{2}$, $C_n=\sqrt {a_n^2+b_m^2}$, y $\theta_n=\arctan \frac {b_n}{a_n}$; y $C_0^*=\frac {a_0^*}{2}$, $C_n^*=\sqrt {a_n^{*2}+b_n^{*2}}$, y $\theta_n=\arctan \frac {b_n^*}{a_n^*}$.

Se aprecia entonces que en cada armónico $n$ , $\theta_n$ determinara el ángulo de desfase en radianes de cada serie de fourier, si queremos obtener el desfase en unidades de tiempo, hay que dividirlo por la frecuencia fundamental ($w_0$), $\frac {\theta_n}{w_0}$, entonces la diferencia $\frac {\theta_n}{w_0} - \frac {\theta_n^*}{w_0}=\frac {\theta_n-\theta_n^*}{w_o}$ determinara el desfase entre los armónicos $n$ de las dos series.

En definitiva, los coeficientes de Fourier también permiten analizar la covarianza cruzada y los desfases que se dan entre las frecuencias relevantes de dos series armónicas.

Ejemplo 5

Se realiza una representacion del cross-espectrum de las series datos anuales del PIB en Indices de Volumen y el consumo de energia final en España, correspondientes al periodo 1995-2018.

Se utiliza la función “crossSpectrum” de la libreria IRISSeismic, la función tiene una opción, spans, donde se puede incluir el vector de enteros impares que se utilizarán para suavizar el periodograma (anchos modificados de Daniels).

En el ejemplo se representan las funciones de espectro de fase y coherencia

library(IRISSeismic)
library(vars)

## Loading required package: MASS

## 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

# Lectura de datos

energia <-read.csv("energia.csv",header=TRUE,sep=";")
E=energia$CEEF
P=energia$PIB_IV

ts1=ts(P,frequency = 1, start = 1995)
ts2=ts(E,frequency = 1, start = 1995)


# Calculate the cross spectrum
DF <- crossSpectrum(ts.union(ts1,ts2),spans=4)

# Calculate the transfer function
transferFunction <- DF$Pxy / DF$Pxx
transferAmp <- Mod(transferFunction)
transferPhase <- pracma::mod(Arg(transferFunction) * 180/pi,360)

# 2 rows
layout(matrix(seq(2)))

# Plot
plot(1/DF$freq,transferAmp,type='l',log='x',
     xlab="Period",
     main="Transfer Function Amplitude")

plot(1/DF$freq,transferPhase,type='l',log='x',
     xlab="Periodo", ylab="degrees",
     main="Transfer Function Phase")

Producto escalar de dos series armónicas.

El producto de dos series armónicas de diferente frecuencia: $x(t)=[a_j\cos (\omega_j)+b_j\sin (\omega_j)]$ e $y(t)=[a_i\cos (\omega_i)+b_i\sin (\omega_i)]$

da como resultado la siguiente suma:

\[\begin{array}{c} a_ja_i\cos(\omega_j)\cos(\omega_i)+a_jb_i\cos (\omega_j)\sin (\omega_i)\\ +a_ib_j\sin (\omega_j)\cos (\omega_i)b_i\sin (\omega_i)+b_jb_i\sin(\omega_j)\sin(\omega_i) \end{array}\]

considerando las identidades del producto de senos y cosenos, quedaría:

\[\begin{array}{c} \frac{a_ja_i+b_jb_i}{2} \cos(\omega_j- \omega_i)+\frac{b_ja_i-b_ja_j}{2}\sin(\omega_j- \omega_i)\\ +\frac{a_ja_i-b_jb_i}{2}\cos(\omega_j+ \omega_i)+\frac{b_ja_i+b_ja_i}{2}\sin(\omega_j+ \omega_i) \end{array}\]

La circularidad de $\omega$ determina que la serie producto de dos series en $t$, resulte una nueva serie cuyos coeficientes de Fourier sean una combinación lineal de los coeficientes de Fourier de las series multiplos.

Partiendo de las dos series siguientes:

\[ \begin{array} {cc} y_t=\eta^y+a_0^y\cos(\omega_0)+b_0^y\sin(\omega_0)+a_1^y\cos(\omega_1)+b_1^y\sin(\omega_1)+ a_2^y\cos(\omega_2)+b_2^y\sin(\omega_2)+a_3^y\cos(\omega_3)\\ x_t=\eta^x+a_0^x\cos(\omega_0)+b_0^x\sin(\omega_0)+a_1^x\cos(\omega_1)+b_1^x\sin(\omega_1)+ a_2^x\cos(\omega_2)+b_2^x\sin(\omega_2)+a_3^x\cos(\omega_3) \end{array} \]

Dada una matriz $\Theta^{\dot x\dot x}$ de tamaño 8x8 :

\[ \Theta^{\dot x\dot x} = \eta^x I_8+\frac{1}2\left( \begin{array}{cccccccc} 0& a_0^x& b_0^x & a_1^x & b_1^x & a_2^x & b_2^x& 2a_3^x \\ 2a_0^x& a_1^x& b_1^x & a_0^x+a_2^x & b_0^x+b_2^x & a_1^x+2a_3^x & b_1^x& 2a_2^x \\ 2b_0^x& b_1^x&- a_1^x & -b_0^x+b_2^x & a_0^x-a_2^x &- b_1^x &a_1^x- a_3^x &- 2b_2^x \\ 2a_1^x& a_0^x+a_2^x&- b_0^x+b_2^x & 2a_3^x &0 & a_0^x+a_2^x & b_0^x-b_2^x& 2a_1^x \\ 2b_1^x& a_0^x+b_2^x&- b_0^x-a_2^x &0& -2a_3^x & -b_0^x+b_2^x & a_0^x-a_2^x& -2b_1^x \\ 2a_2^x& a_1^x+2a_3^x&- b_1^x & a_0^x+a_2^x &-b_0^x-b_2^x & a_1^x &- b_1^x& 2a_0^x \\ 2b_2^x& b_1^x& a_1^x-2a_3^x & b_0^x-b_2^x &a_0^x-a_2^x & -b_1^x &- a_1^x& -2b_0^x \\ 2a_3^x& a_2^x& -b_2^x & a_1^x &- b_1^x & a_0^x & -b_0^x& 0 \end{array} \right) \] Se demuestra que:

$\dot z=\Theta^{\dot x\dot x}\dot y$

donde $\dot y = Wy$,$\dot x = Wx$, y $\dot z = Wz$.

En el dominio del tiempo:

$z_t= x_t y_t=W^T\dot x W^T\dot y=W^T Wx_t W^T\dot y=x_tI_nW^T\dot y$

$W^T\dot z=x_tI_nW^T\dot y$

$\dot z=Wx_tI_nW^T\dot y$

Entonces:

$\Theta^{\dot x\dot x}=W^Tx_tI_nW$

La matriz cuadrada $\Theta^{\dot x\dot x}$ puede ser utilizada para obtener los resultados en el dominio de la frecuencia de diversas funciones de series de tiempo . Por ejemplo, si se desea obtener el desarrollo de los coeficientes en fourier de $z_t=x_t^2$, entonces:

$\dot z= Wx_tI_nW^T\dot x$

En consecuencia, si $z_t=x_t^n$

$\dot z= Wx_t^{n-1}I_nW^T\dot x$

Si ahora queremos obtener el desarrollo en coeficientes de fourier de $z_t=\frac{x_t}{y_t}$, entonces:

$\dot z= W[\frac{1}y_t]I_nW^T\dot x$

Función cdf(a)

Obtiene la matriz auxiliar para operaciones con vectores en dominio de tiempo y dominio de la frecuencia, pre-multiplica un vector por la matriz ortogonal, $W$ y por su transpuesta, Parra F. (2013)

cdf <- function(y) {
# Author: Francisco Parra Rodríguez 
# http://econometria.wordpress.com/2013/08/21/estimation-of-time-varying-regression-coefficients/ 
a <- matrix(y, nrow=1)
n <- length(y)
uno <- as.numeric (1:n)
A <- MW(n)
I<- diag(c(a))
B <- A%*%I
B%*%t(A)
}

Aproximación bivariada con series de Fourier: Forma Flexible de Fourier (FFF)

Gallant (1981;1982) introdujo una forma funcional con capacidades muy distintas a las propuestas hasta el momento; cuyas propiedades de flexibilidad eran en todos los casos locales. La forma de Fourier que utiliza Gallant posee la propiedad de flexibilidad global; es decir; 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.

Por tratarse de una forma Sobolev-flexible (frente a la Diewert-flexibilidad de las anteriores) es capaz de estimar consistentemente las elasticidades precio y renta sobre todo el espacio de datos (ElBadawi, Gallant y Souza; 1983); además; asintóticamente pueden conseguirse contrastes estadísticos insesgados (Gallant; 1981 y 1982) y la eliminación del problema de inferencias aumentadas provocado por la especificación de un determinado modelo. Por último; Gallant y Souza (1991) han mostrado la normalidad asintótica de las estimaciones derivadas de la forma de Fourier.

En la parte negativa, el modelo de Fourier puede conseguir la regularidad global, pero las restricciones paramétricas que ello implica son excesivamente fuertes (Gallant, 1981); sin embargo, existen condiciones más débiles (que no destruyen ni la flexibilidad ni la consistencia de los estimadores) con las que se puede conseguir la regularidad teórica al menos sobre un conjunto finito de puntos (Gallant y Golub, 1983); aunque la implementación de tales restricciones resulta compleja (McFadden; 1985). En cualquier caso, las simulaciones de Monte Carlo realizadas por Fleissig, Kastens y Terrell (1997) y Chalfant y Gallant (1985) han mostrado que la región de regularidad de la forma de Fourier libre -sin restricciones de ningún tipo- es mucho mayor que la correspondiente a las formas Leontief-Generalizada o Translog.

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

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

La aproximación a una función no periódica g(x) por una serie de expansión de Fourier se realiza en Gallant (1981) añadiendo es esta un término lineal y cuadrático. De esta forma que la aproximación univariada se escribe como:

\[g(x/\theta)=\alpha + \beta x + \frac 1 2 \delta x^2 + 2\sum_{j=1}^J[a_j\cos(jx)+b_j\sin(jx)] (3)\]

El vector de parámetros es $(\theta)=(\alpha,\beta,\delta ,a_1,...,a_J, b_1,....,b_J)$ de longitud $K=3+2J$, siendo siendo $J \approx \sqrt n$ .

La expresión de la primera y segunda derivada de la función (3) son las siguientes:

\[D_xg(x/\theta)=\beta + \delta x + \sum_{j=1}^J[-a_j\sin(jx)-b_j\cos(jx)]j\]

\[D_x^2g(x/\theta)=\delta + \sum_{j=1}^J[-a_j\cos(jx)+b_j\sin(jx)]j^2\]

Dado que la variable exógena $x$ no está expresada en forma periódica, debe de transformase o normalizarse en un intervalo de longitud menor que $2\pi$;$[0,2\pi ]$.

Función FFF(y,x)

La función realiza la aproximación de la expansión en series de fourier descrita en Gallant (1981) para una variable dependiente (y) y otra independiente (x).

# Funcion FFF
FFF=function(y,x){
z=2*pi*x/max(x)
n=3+2*abs(sqrt(length(x)))
X=data.frame(X=x,X2=(x^2)/2,X3=cos(z),X4=sin(z))
m=trunc((n-5)/2)
for(i  in (1:m))
{X5=cos(2*i*z)
X6=sin(2*i*z)
X=cbind(X,X5,X6)}
list(fitted=lm(y~as.matrix(X))$fitted,X=X,residuals=lm(y~as.matrix(X))$residuals)}

Ejemplo 6

A continuación va a realizarse una FFF con datos anuales del PIB en Indices de Volumen y el consumo de energia final en España, correspondientes al periodo 1995-2018.

Dado que la longitud de los datos del ejemplo son 24, la expansión considera que $K=3+2*5=13$

energia <-read.csv("energia.csv",header=TRUE,sep=";")
E=energia$CEEF
P=energia$PIB_IV
fff.E=FFF(E,P)
plot(ts(E,frequency = 1, start = 1995),type="l",main="Consumo Energia Final (Ktep).España",ylab="")
lines(ts(lm(E~P)$fitted, frequency = 1, start = 1995), type="l", col=2)
lines(ts(fff.E$fitted, frequency = 1, start = 1995), type="l", col=3)
legend("top", ncol=3, c("CEF","Estimado LM","Estimado FFF"), cex=0.6, bty="n", fill=c(1,2,3))

Aproximación FFF multivariada:

La aproximación multivariada se describe en Gallant (1984):

\[g(x/\theta)=u_0 + b' x + \frac 1 2 x'Cx + \sum_{\alpha=1}^A\left[u_{0\alpha}+2 \sum_{j=1}^J m_{j\alpha}\cos(jk'_{\alpha}z)+n_{j\alpha}\sin(jk'_{\alpha}z)\right]\]

Donde, $x$ es un vector de $Nx1$ variables, $b$ es un vector de $Nx1$ coeficientes, $C$ es una matriz simétrica de de $NxN$ coeficientes, $z$ es un vector $Nx1$ de valores transformados de $x$ ; $m_{j\alpha}$ y $n_{j\alpha}$ son coeficientes, $k'_{\alpha}=[k_{x_1},k_{x_2},...,k_{x_N}]$ son multi-índices; vectores de $1xN$ elementos que representan a la derivadas parciales de una función de producción de Fourier para los diferentes tipos de expansión.

Dado que las variable exógenas no están expresada en forma periódica, deben de transformase o normalizarse en un intervalo de longitud menor que $2\pi$. Fulginiti et all (2003) sugieren transformar el vector de variables $x$ definiendo:

\[l_i =log(a_i)+log(x_i) ; i = 1,2,..., N\] siendo $a_i=Min[log(x_i)]+10^{-5}$

Entonces el valor de las variable transformada sería:

\[ z_{it}=j\lambda k'_{\alpha}[log(a_t)+log(x_t)]\]

siendo $\lambda=\frac {2\pi - \epsilon}{max(l_i)}$

donde $\epsilon$ es valor positivo arbitrario y pequeño, si bien se recomienda escoger: $2\pi - \epsilon=6$.

La forma flexible que aproxima a una ecuación con tres variables exógenas, y $j=1$ sería:

\[g(x/\theta)=u_0 + \sum_{i=1}^3 b_ix_{it} + \frac 1 2 \sum_{i=1}^3 c_{ii}x_{it}^2 + \sum_{i=1}^3 \sum_{k>1}^3 c_{ik}x_{it} x_{kt}+ [m_{t} cos(z_t) + n_{t} sin (z_t)]+\] \[+\sum_{i=1}^3 \sum_{k>1}^3 [m_{ik} cos (z_{it}-z_{kt})+ n_{ik} sin (z_{it}-z_{kt})]+\sum_{i=1}^3 \sum_{k>1}^3 [m_{ik1} cos (z_{it}-z_{kt}-z_t)+n_{ik1} sin (z_{it}-z_{kt}-z_t)]+\] \[+\sum_{i=1}^3 \sum_{k>1}^3 [m_{ik2} cos (z_{it}-z_{kt}+z_t)+n_{ik2} sin (z_{it}-z_{kt}+z_t)]\]

Es habitual en este tipo de trabajos incluir una tendencia temporal como variable exógena, en cuyo caso la forma flexible sería:

\[g(x/\theta)=u_0 + \sum_{i=1}^3 b_ix_{it} + \frac 1 2 \sum_{i=1}^3 c_{ii}x_{it}^2 + \sum_{i=1}^3 \sum_{k>1}^3 c_{ik}x_{it} x_{kt}+ b_t + \frac 1 2 b_{tt} t^2 +\sum_{i=1}^3 b_{it}x_{it}t + [m_{t} cos(z_t) + n_{t} sin (z_t)]+ \] \[+\sum_{i=1}^3 \sum_{k>1}^3 [m_{ik} cos (z_{it}-z_{kt})+ n_{ik} sin (z_{it}-z_{kt})]+\sum_{i=1}^3 \sum_{k>1}^3 [m_{ik1} cos (z_{it}-z_{kt}-z_t)+n_{ik1} sin (z_{it}-z_{kt}-z_t)]+\]

\[+\sum_{i=1}^3 \sum_{k>1}^3 [m_{ik2} cos (z_{it}-z_{kt}+z_t)+n_{ik2} sin (z_{it}-z_{kt}+z_t)]\] ## Ejemplo 7

La función de Coob-Douglas una de las funciones de producción más empleadas en el ámbito de la economía, la función fue estimada por Paul Douglas y Charles Cobb en 1927. Se expresa:

\[y_t=Ak_t^{\alpha}l_t^{\beta}\] siendo $y_t$ la producción nacional, $k_t$ el stock de capital, y $l_t$ el empleo.

Con los datos de la economía americana de 1899 a 1922, estimamos la aproximación FFF multivariada:

\[Ln(y_t)=\mu_o + b_1 Ln(k_t) + b_2 Ln(l_t) + b_{11} (Ln(k_t))^2 + b_{22} (Ln(l_t))^2 + b_{12}(Ln(k_t)Ln(l_t))+ [m1 cos(k*_t) + n1 sin(k*_t)] + [m2 cos(l*_t) + n2 sin(l*_t)] + [m_{11} cos(2k*_t) + n_{11} sin(2k*_t)] + [m_{22} cos(2l*_t) + n_{22} sin(2l*_t)] + [m_{12} cos(k*_t+l*_t) + n_{12} sin(k*_t + l*_t)] \]

datos <-read.csv("coob_douglas.csv",header=TRUE,sep=";")
y=log(datos$Y)
k=log(datos$K)
l=log(datos$L)
k.2=k^2
l.2=l^2
kl.2=k*l
k2=2*pi*k/max(k)
l2=2*pi*l/max(l)
k.3.1=cos(k2)
k.3.2=sin(k2)
l.3.1=cos(l2)
l.3.2=sin(l2)
k.4.1=cos(2*k2)
k.4.2=sin(2*k2)
l.4.1=cos(2*l2)
l.4.2=sin(2*l2)
kl.4.1=cos(k2+l2)
kl.4.2=sin(k2+l2)
X=data.frame(y,k,l,k.2,l.2,kl.2,k.3.1,k.3.2,l.3.1,l.3.2,k.4.1,k.4.2,l.4.1,l.4.2,kl.4.1,kl.4.2)
mod1=lm(y~.,data=X)
summary(mod1)

## 
## Call:
## lm(formula = y ~ ., data = X)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.037278 -0.004165  0.000928  0.009010  0.020974 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.208e+05  2.463e+05   1.709  0.12588   
## k           -1.735e+04  3.969e+03  -4.370  0.00238 **
## l           -1.687e+05  1.082e+05  -1.560  0.15740   
## k.2          1.593e+03  3.662e+02   4.351  0.00244 **
## l.2          1.737e+04  1.091e+04   1.593  0.14986   
## kl.2         2.634e+01  1.322e+01   1.993  0.08143 . 
## k.3.1        3.119e+03  6.972e+02   4.473  0.00207 **
## k.3.2       -2.569e+03  6.163e+02  -4.168  0.00313 **
## l.3.1        2.937e+04  1.923e+04   1.528  0.16514   
## l.3.2       -1.703e+04  8.298e+03  -2.053  0.07421 . 
## k.4.1       -5.488e+01  1.403e+01  -3.912  0.00447 **
## k.4.2        2.675e+02  6.249e+01   4.281  0.00268 **
## l.4.1       -1.189e+03  9.406e+02  -1.265  0.24163   
## l.4.2        1.983e+03  9.824e+02   2.019  0.07819 . 
## kl.4.1       1.711e+01  7.908e+00   2.163  0.06250 . 
## kl.4.2      -1.332e+01  9.337e+00  -1.427  0.19155   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02246 on 8 degrees of freedom
## Multiple R-squared:  0.9976, Adjusted R-squared:  0.993 
## F-statistic: 219.8 on 15 and 8 DF,  p-value: 8.989e-09

plot(ts(exp(y),frequency = 1, start = 1889),type="l",main="Y",ylab="")
lines(ts(exp(lm(y~k+l)$fitted), frequency = 1, start = 1889), type="l", col=2)
lines(ts(exp(mod1$fitted), frequency = 1, start = 1889), type="l", col=3)
legend("top", ncol=3, c("Y","Estimado LM","Estimado FFF"), cex=0.6, bty="n", fill=c(1,2,3))

Aproximación multivariada FFF utilizando funciones paramétricas

La aproximación FFF multivariada de Gallant (1981,1983) presenta dificultades prácticas ya que precisa de una gran cantidad de datos para ser estimada por los métodos convencionales, la reducción de grados de libertad que ocasiona el utilizar secuencias de series de senos y cosenos en una regla de difícil aplicación práctica, puede solventarse parametrizando los ángulos que determinan la relación polar en un eje de tres dimensiones.

Se denominan ecuaciones paramétrica a aquellas ecuaciones en que las variables $X$ e $Y$, cada una separadamente, están expresadas en función de la misma tercera variable, $t$, a la que se denomina variable paramétrica, estas ecuaciones se representan en la siguiente forma general:

\[\begin{cases} & \ X=u(t) \\ & \ Y=v(t) \end{cases}\]

Una ecuación paramétrica permite representar curvas o superficies en el plano o en el espacio, mediante valores arbitrarios (parámetros). En el uso estándar del sistema de coordenadas, una o dos variables (dependiendo de si se utilizan dos o tres dimensiones respectivamente) son consideradas como variables independientes, mientras que la restante, a la que se denomina variable dependiente, toma un valor en función de los valores que toman las variable(s) independiente(s). Así por ejemplo la expresión de un punto cualquiera $(x,y)$ equivale a la expresión $(x,f (x))$.

Esta representación tiene la limitación de requerir que la curva sea una función de $X$ en $Y$, es decir que todos los valores X tengan un valor y sólo un valor correspondiente en $Y$, y no todas las curvas cumplen con dicha condición. En una ecuación paramétrica, tanto $X$ como $Y$ son considerados variables dependientes, cuyo resultado surge de una tercera variable (sin representación gráfica) conocida como parámetro, lo que la representación de funciones circulares en donde un valor de $X$ puede dar lugar a dos valores de $Y$.

Por ejemplo, la posible parametrización de la expresión $Y=X^2$, de la forma $\begin{cases} & \ X=u(t) \\ & \ Y=v(t) \end{cases}$, sería $\begin{cases} & \ X=t \\ & \ Y=t^2 \end{cases}$.

Una circunferencia con centro en el origen de coordenadas y radio $r$ verifica que $X^2+Y^2=r^2$, y su expresión paramétrica sería $\begin{cases} & \ X=rcos(t) \\ & \ Y=rsin(t) \end{cases}$.

La representación paramétrica de una curva en un espacio n-dimensional consiste por tanto en “n” funciones de una variable t que actúa como variable independiente o parámetro (habitualmente se considera que t es un número real y que los puntos del espacio n-dimensional están representados por n coordenadas reales), de la forma:

\[ e_i=f_i(t),f_i:[a.b] \rightarrow \Re \] donde $e_i$ representa la i-ésima coordenada del punto generado al asignar valores del intervalo $[a, b]$ a $t$. Por ejemplo, para representar una curva en el espacio se usan 3 funciones $X = u(t)$ , $Y = v(t)$ y $Z = g(t)$.

Es común que se exija que el intervalo $[a, b]$ sea tal que a cada punto $a \leq y \leq b$ le corresponda un punto distinto de la curva; si las coordenadas del punto obtenido al hacer $t = a$ son las mismas del punto correspondiente a $t =b$ la curva se denomina cerrada.

Se dice que un punto de la curva correspondiente a un valor t del intervalo es un punto ordinario si las derivadas de las funciones paramétricas existen en y son continuas en ese punto y al menos una es distinta de 0. Si un arco de curva está compuesto solamente de puntos ordinarios se denomina suave. Es común resumir las ecuaciones paramétricas de una curva en una sola ecuación vectorial:

\[\vec r_t=\sum_{i=1}^n f_i(t)\hat e_i=f_1(t)\hat e_1+f_2(t) \hat e_2 +...+f_n(t)\hat e_n\] donde $\hat e_i$ representa al vector unitario correspondiente a la coordenada i-ésima.

Por ejemplo, las funciones paramétricas de un círculo unitario con centro en el origen son $X = cos (t)$, $Y = sen (t)$. Podemos reunir estas ecuaciones como una sola ecuación de la forma:

\[\vec r_t=cos(t)\hat i +sin(t) \hat j \] Una superficie parametrizada en $\Re^3$ es la imagen de una función continua $S$ definida en una región $D \subseteq \Re ^2$ que toma valores en $\Re^3$, esto es,

\[S:(u,v) \in D \subseteq \Re ^2 \rightarrow S(n,v) = (x(u,v), y(u,v), z(u,v)) \subseteq \Re^3\]

Las variables independientes de la función $S$ se llaman parámetros de la superficie y la propia función $S$ recibe el nombre de parametrización de la superficie. La imagen por $S$ de la frontera de la región $D$ se llama borde o contorno de la superficie. Si $S$ es inyectiva, lo que significa que no hay puntos dobles, entonces se dice que la superficie es simple.

Dadas $n$ observaciones de la variable aleatoria $x_t \rightarrow N(\mu_x,\sigma_x^2)$, a cada obseración $(t,x_t)$ le corresponde un punto $P_t$ en el eje cartesiano, de forma que:

\[P_1 \Rightarrow (1,x_1), P_2 \Rightarrow (2,x_2),...,P_n \Rightarrow (n,x_n) \] que a su vez, se hace corresponder con una forma polar $r_t$, para cada par $(t,x_t)$, siendo:

$r_t=\sqrt {t^2+x_t^2}$ y $\gamma_t=arctan \frac {x_t}{t}$

Dado $P_t(t,x_t)=t+ix_t=[r_tcos(\gamma_t)]+i[r_tsin(\gamma_t)]$, se deduce que $x_t=r_tsin(\gamma_t)$ y $x_t=tan(\gamma_t)t$.

Si se estima la variable $\gamma_t$ a partir de la expansión FFF de la forma:

\[\gamma_t(t/\theta)=\alpha +\beta t + \frac 1 2 \delta_t^2 + \sum_{j=1}^J [a_jcos(jw_0t) + b_jsin(jw_ot)])\]

siendo $w_0=\frac {2\pi}{n}$, la función $x_t$ quedaría parametrizada en función de $t$ con la siguiente expresión:

\[x_t=tan(\gamma_t)t\]

Tenemos ahora $n$ observaciones de la variable aleatoria $x_t \rightarrow N(\mu_x,\sigma_x^2)$ e $y_t \rightarrow N(\mu_y,\sigma_y^2)$, a cada obseración $(x_t,y_t)$ le corresponde un punto $P_t$ en el eje cartesiano, de forma que:

\[P_1 \Rightarrow (x_1,y_1), P_2 \Rightarrow (x_2,y_2),...,P_n \Rightarrow (x_n,y_n) \] que a su vez, se hace corresponder con una forma polar $r_t$, para cada par $(x_t,y_t)$, siendo:

$r_t=\sqrt {x_t^2+y_t^2}$ y $\alpha_t=arctan \frac {y_t}{x_t}$

Dado $P_t(x_t,y_t)=x_t+iy_t=[r_tcos(\alpha_t)]+i[r_tsin(\alpha_t)]$, se deduce que $y_t=r_tcos(\alpha_t)$ , $x_t=r_tcos(\alpha_t)$ y $y_t=tan(\alpha_t)x_t$.

La variable aleatoria $y_t$ puede parametrizarse como:

\[y_t=tan(\alpha_t)tan(\gamma_t)t\] Las ecuaciones paramétricas serían entonces:

\[\begin{cases} & \ X=x_t=tan(\gamma_t)t \\ & \ Y=y_t=tan(\alpha_t)tan(\gamma_t)t \end{cases}\]

estimandose el ángulo $\alpha_t$ con una expansión FFF, al igual que $\gamma_t$:

\[\alpha_t(t/\theta)=\alpha' +\beta' t + \frac 1 2 \delta_t^{'2} + \sum_{j=1}^J [a'_jcos(jw_0t) + b'_jsin(jw_ot)])\]

Tenemos ahora $n$ observaciones de la variable aleatoria $x_t \rightarrow N(\mu_x,\sigma_x^2)$ , $y_t \rightarrow N(\mu_y,\sigma_y^2)$ y $z_t \rightarrow N(\mu_z,\sigma_z^2)$. En este caso las relaciones geométricas a considerar son las que aparecen en la figura adjunta:

Se parte ahora de la representación polar entre cada par $x_t$ y $z_t$, que vendrá dada por el modulo $r_t=\sqrt {x_t^2+z_t^2}$ y el argmento $\alpha_t=arctan \frac {z_t}{x_t}$. Construimos ahora un nuevo plano entre $r_t$ e $y_t$.

Dado que el modulo $r_t$ puede tener un valor diferente según se cambie el nivel de la variable parece aconsejable normalizar las variables.

En consecuencia ahora tenemos dos variables $r_t$ e $y_t$ cuya representación polar tendrá a su vez un módulo $\varphi_t=\sqrt {r_t^2+y_t^2}$ y el argumento $\beta_t=arctan \frac {y_t}{r_t}$.

Operando $\varphi_t=\sqrt {z_t^2+x_t^2+y_t^2}$, y la representación polar del sistema vendría dada por : $r_t=\varphi_t cos (\beta_t)$,$y_t=\varphi_t sin (\beta_t)$, e $y_t=tan (\beta_t) \sqrt {x_t^2+z_t^2}$.

Dado que $z_t=tan (\alpha_t)x_t$, entonces:

\[y_t=tan (\beta_t) \sqrt {x_t^2+[tan (\alpha_t)x_t]^2}=x_t\sqrt {tan^2 (\beta_t)[1+tan^2 (\alpha_t)]}\] considerando tanto la sucesión de ángulos $\beta_t$ y $\alpha_t$ como series de Fourier, se puede afirmar que el conjunto de datos $(x_t,y_t,z_t)$, puede parametrizarse en función de una de ellas cualesquiera, y $t$.

Supongamos que en nuestro conjunto de datos la dimensión $x_t$, es exógena, entonces:

\[\begin{cases} & \ X=x_t \\ & \ Z=z_t=tan(\alpha_t)x_t \\ & \ Y=y_t=\sqrt {tan^2 (\beta_t)[1+tan^2 (\alpha_t)]} x_t \end{cases}\]

siendo:

\[\beta_t(t/\theta)=\alpha^* +\beta^* t + \frac 1 2 \delta_t^{*2} + \sum_{j=1}^J [a_j^*cos(jw_0t) + b_j^*sin(jw_ot)])\]

La parametrización del sistema a la dimensión $z_t$:

\[\begin{cases} & \ Z=z_t \\ & \ X=x_t=tan(\frac {\pi}{2} - \alpha_t)z_t \\ & \ Y=y_t=\sqrt {tan^2 (\beta_t)[1+tan^2 (\frac {\pi}{2} - \alpha_t)x_t]} z_t \end{cases}\]

ya que $\frac {\pi}{2} - \alpha_t=arctan \frac {x_t}{z_t}$.

Por ultimo la parametrización sobre $y_t$:

\[\begin{cases} & \ Y=y_t \\ & \ X=x_t=\frac {y_t}{\sqrt {tan^2 (\beta_t)[1+tan^2 (\alpha_t)]}} \\ & \ Z=z_t=\frac {y_t}{\sqrt {tan^2 (\beta_t)[1+tan^2 (\frac {\pi}{2} - \alpha_t)]}} \end{cases}\]

Considerando la parametrización en $t$ de $x_t$, el sistema quedaría:

\[\begin{cases} & \ X=x_t=tan(\gamma_t)t \\ & \ Z=z_t=tan(\alpha_t)tan(\gamma_t)t \\ & \ Y=y_t=\sqrt {tan^2 (\gamma_t)tan^2 (\beta_t)[1+tan^2 (\alpha_t)} t \end{cases}\]

Ejemplo 8

Utilizando los datos del Ejemplo 7, estimamos una función FFF para el argumento $\gamma_t$.

# Calculo del angular
t=seq(1:length(y))
g=atan(l/t)
# Expansión del angular
fff.g=FFF(g,t)
# Estimación de L
l.fitted=tan(fff.g$fitted)*t
# Representación gráfica
plot(ts(g,frequency=1,start=1899),type="l",main="Angular",ylab="")
lines(ts(fff.g$fitted,frequency=1,start=1899),type="p",col=2)

plot(ts(exp(l),frequency=1,start=1899),type="l",main="Empleo",ylab="")
lines(ts(exp(l.fitted),frequency=1,start=1899),type="p",col=2)

La estimación del stock de capital a partir del empleo, utilizando funciones paramétricas y una expansión FFF para el argumento $\alpha_t$:

# Calculo del angular
a=atan(k/l)
# Expansión del angular
fff.a=FFF(a,t)
# Estimación de K
k.fitted=tan(fff.a$fitted)*l
# Representación gráfica
plot(ts(a,frequency=1,start=1899),type="l",main="Angular",ylab="")
lines(ts(fff.a$fitted,frequency=1,start=1899),type="p",col=2)

plot(ts(exp(k),frequency=1,start=1899),type="l",main="Stock Capital",ylab="")
lines(ts(exp(k.fitted),frequency=1,start=1899),type="p",col=2)

# Calculo del angular
r=sqrt(l^2+k^2)
b=atan(y/r)
# Expansión del angular
fff.b=FFF(b,t)
# Estimación de Y
y.fitted=sqrt(1+(tan(fff.a$fitted))^2)*l*tan(fff.b$fitted)
# Representación gráfica
plot(ts(b,frequency=1,start=1899),type="l",main="Angular",ylab="")
lines(ts(fff.b$fitted,frequency=1,start=1899),type="p",col=2)

plot(ts(exp(y),frequency=1,start=1899),type="l",main="Producción",ylab="")
lines(ts(exp(y.fitted),frequency=1,start=1899),type="p",col=2)

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.

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 \ \ \ (1)\]

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 (1) del dominio del tiempo al dominio de la frecuencia en Engle (1974), a partir de la matriz $W$ (Ver Análisis Espectral), 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 $x$ $y$ por $W$(ver apartado), obtenemos:

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

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

Si el vector de las perturbaciones en (1) 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 (2) 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 (1).

Estimar (2) 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.

Ejemplo 9

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)

## 
## Attaching package: 'descomponer'

## The following objects are masked _by_ '.GlobalEnv':
## 
##     cdf, gdf, gdt, gperiodograma, gtd, MW, periodograma, td

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 \ \ \ \ \ \ (3)\]

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

Premultiplicando por $W$ obtenemos:

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

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

El sistema (4) puede reescribirse como:

\[ \dot y=Wx_tI_nW^T\dot \beta + WI_nW^T\dot u\]

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.

Ejemplo 10

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:

Convertimos al dominio de tiempo cada regresor (cada columna de la matriz X)
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
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
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
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
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
y sucesivamente ….
Se combinan las series pronosticadas con los coeficientes que resultan de la regresión en el dominio de la frecuencia.

Ejemplo 11

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.

Funcion MW2 (n,m)

La función MW2, obtiene la Matriz $W$ desfasada m periodos de una serie temporal de tamaño n.

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 gdf2(y,m)

La función gdf2, transforma al dominio de la frecuencia la serie temporal $y$ adelanta en $m$ periodos.

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

función gdt2(y,m)

La función gdf2, transforma al dominio del tiempo la serie temporal $y$ presentada en el dominio de la frecuencia, y adelantada en $m$ periodos.

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

Ejemplo 12

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, la función “rdf” permite realizar dicha regresión.

funcion rdf(x,y)

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

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}\]

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.
Calcula la matriz $Wx_tI_nW^T$ y la ordena en base al indice anterior.
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.
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$.

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

rdf <- function (y,x) {
  # Author: Francisco Parra Rodriguez
  # http://rpubs.com/PacoParra/24432
  # Leemos datos en forma matriz
  a <- matrix(y, nrow=1)
  b <- matrix(x, nrow=1)
  n <- length(a)
  # calculamos el cros espectro mediante la funcion cperiodograma
  cperiodograma <- function(y,x) {
    # Author: Francisco Parra Rodriguez
    # http://econometria.wordpress.com/2013/08/21/estimation-of-time-varying-regression-coefficients/ 
    cfx <- gdf(y)
    n <- length(y)
    cfy <- gdf(x)
    if (n%%2==0) {
      m1x <- c(0)
      m2x <- c()
      for(i in 1:n){
        if(i%%2==0) m1x <-c(m1x,cfx[i]) else m2x <-c(m2x,cfx[i])}
      m2x <- c(m2x,0)
      m1y <- c(0)
      m2y <- c()
      for(i in 1:n){
        if(i%%2==0) m1y <-c(m1y,cfy[i]) else m2y <-c(m2y,cfy[i])}
      m2y <-c(m2y,0) 
      frecuencia <- seq(0:(n/2)) 
      frecuencia <- frecuencia-1
      omega <- pi*frecuencia/(n/2)
      periodos <- n/frecuencia
      densidad <- (m1x*m1y+m2x*m2y)/(4*pi)
      tabla <- data.frame(omega,frecuencia, periodos,densidad)
      tabla$densidad[(n/2+1)] <- 4*tabla$densidad[(n/2+1)]
      data.frame(tabla[2:(n/2+1),])}
    else {m1x <- c(0)
    m2x <- c()
    for(i in 1:(n-1)){
      if(i%%2==0) m1x <-c(m1x,cfx[i]) else m2x <-c(m2x,cfx[i])}
    m2x <-c(m2x,cfx[n])
    m1y <- c(0)
    m2y <- c()  
    for(i in 1:(n-1)){
      if(i%%2==0) m1y <-c(m1y,cfy[i]) else m2y <-c(m2y,cfy[i])}
    m2y <-c(m2y,cfy[n])
    frecuencia <- seq(0:((n-1)/2)) 
    frecuencia <- frecuencia-1
    omega <- pi*frecuencia/(n/2)
    periodos <- n/frecuencia
    densidad <- (m1x*m1y+m2x*m2y)/(4*pi)
    tabla <- data.frame(omega,frecuencia, periodos,densidad)
    data.frame(tabla[2:((n+1)/2),])}
  }
  cper <- cperiodograma(a,b)
  S1 <- data.frame(f1=cper$frecuencia,p=abs(cper$densidad))
  S <- S1[with(S1, order(-p)), ]
  id <- seq(2,n)
  fpart = function(vec) {
    ret = vec - as.integer(vec)
    ret
  }
  evens=function(vec) {
    stopifnot(class(vec)=="integer")
    ret = vec[fpart(vec/2)==0]
    ret
  }
  odds=function(vec) {
    stopifnot(class(vec)=="integer")
    ret = vec[fpart(vec/2)!=0]
    ret
  }
  m1 <- cbind(S$f1*2,evens(id))
  if (n%%2==0) {m2 <- cbind(S$f1[1:(n/2-1)]*2+1,odds(id))} else 
  {m2 <- cbind(S$f1*2+1,odds(id))}
  m <- rbind(m1,m2)
  colnames(m) <- c("f1","id")
  m <- data.frame(m)
  M=m[with(m, order(id)), ]
  M <- rbind(c(1,1),M)
  # Obtenemos la funcion auxiliar (cdf) del predictor y se ordena segun el indice de las mayores densidades absolutas del co-espectro.
  cx <- cdf(b)
  id <- seq(1,n)
  S1 <- data.frame(cx,c=id)
  S2 <- merge(M,S1,by.x="id",by.y="c")
  S3=S2[with(S2, order(f1)), ]
  m <- n+2
  X1 <- S3[,3:m]
  X1 <- rbind(C=c(1,rep(0,(n-1))),S3[,3:m])
  # Se realizan las regresiones en el dominio de la frecuencia utilizando un modelo con constante, pendiente y los armonicos correspondientes a las frecuencias mas altas de la densidad del coespectro. Se realiza un test de durbin para el residuo y se seleccionan aquellas que son significativas. 
  par <- evens(id)
  i <- 1
  D <- 1
  resultado <- cbind(i,D)
  for (i in par) {
    X <- as.matrix(X1[1:i,])
    cy <- gdf(a)
    B1 <- solve(X%*%t(X))%*%(X%*%cy)
    Y <- t(X)%*%B1 
    F <- gdt(Y)
    res <- (t(a) - F)
    T <- td(res)
    L <- as.numeric(c(T$min<T$s2,T$s2<T$max))
    LT <- sum(L)
    D <- LT-(n-1)
    resultado1 <- cbind(i,D)
    resultado <- rbind(resultado,resultado1)
    resultado}
  resultado2 <-data.frame(resultado)
  criterio <- resultado2[which(resultado2$D==0),]
  sol <- as.numeric(is.na(criterio$i[1]))
  if (sol==1) {
    X <- as.matrix(X1[1:2,])
    cy <- gdf(a)
    B1 <- solve(X%*%t(X))%*%(X%*%cy)
    Y <- t(X)%*%B1 
    F <- gdt(Y)
    res <- (t(a) - F)
    datos <- data.frame(cbind(t(a),t(b),F,res))
    colnames(datos) <- c("Y","X","F","res")
    list(datos=datos,Fregresores=t(X),Tregresores= t(MW(n))%*%t(X),Nregresores=criterio$i[1])  
  } else {
    X <- as.matrix(X1[1:criterio$i[1],])
    cy <- gdf(a)
    B1 <- solve(X%*%t(X))%*%(X%*%cy)
    Y <- t(X)%*%B1 
    F <- gdt(Y)
    res <- (t(a) - F)
    datos <- data.frame(cbind(t(a),t(b),F,res))
    colnames(datos) <- c("Y","X","F","res")
    sse=sum(res^2)
    gcv=length(res)*sse/((length(res)-dim(X)[1])^2)
    list(datos=datos,Fregresores=t(X),Tregresores= t(MW(n))%*%t(X),Nregresores=criterio$i[1],sse=sse,gcv=gcv)}
}

funcion rdf2(y,x)

Realiza una aproximación entre la serie x e y con la tecnica RBS, utilizando como regresores una tendencia cuadráticas y los armónicos correspondientes a las $K=\frac {\sqrt n} 2$ frecuencias mas bajas.

Dado que la variable exógena $x$ no está expresada en forma periódica, debe de transformase o normalizarse en un intervalo de longitud menor que $2\pi$;$[0,2\pi ]$.

rdf2 <- function (y,x) {
  # Author: Francisco Parra Rodríguez
    # Leemos datos en forma matriz
    a <- matrix(y, nrow=1)
    b <- matrix(x, nrow=1)
  n=ifelse(length(a)%%2==0,round(2*sqrt(length(a))),round(2*sqrt(length(a)-1)))
    # Obtenemos la funcion auxiliar (cdf) del predictor y se ordena segun el indice de las mayores densidades absolutas del co-espectro.
  cx <- cdf(b)
  cx <- cx[1:(n+1),]
  cx.2=t(gdf(b^2))/10
  X <- rbind(C=c(1,rep(0,(n-1))),cx,cx.2)
  # Se realizan las regresiones en el dominio de la frecuencia utilizando un modelo con constante, pendiente y los arm?nicos correspondientes a las raiz de n frecuencias mas altas
  cy <- gdf(a)
  B1 <- solve(X%*%t(X))%*%(X%*%cy)
  Y <- t(X)%*%B1 
  F <- gdt(Y)
  res <- (t(a) - F)
  datos <- data.frame(cbind(t(a),t(b),F,res))
  colnames(datos) <- c("Y","X","F","res")
list(datos=datos,Fregresores=t(X),Tregresores= t(MW(length(a)))%*%t(X),Nregresores=n+3,Betas=B1)}

Ejemplo 13

A continuación va a realizarse una RBS con con datos anuales del PIB en Indices de Volumen y el consumo de energia final en España, correspondientes al periodo 1995-2018. 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).

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

## 
## Call:
## lm(formula = E ~ P)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2040.3  -415.2   117.2   503.3  1314.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2395.15    1266.05  -1.892   0.0718 .  
## P             229.49      13.71  16.741 5.28e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 827.2 on 22 degrees of freedom
## Multiple R-squared:  0.9272, Adjusted R-squared:  0.9239 
## F-statistic: 280.3 on 1 and 22 DF,  p-value: 5.285e-14

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

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

## 
## Call:
## lm(formula = E.1 ~ 0 + K.1 + P.1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3143.99   -37.02   149.20   307.14  1519.82 
## 
## Coefficients:
##     Estimate Std. Error t value Pr(>|t|)    
## K.1 -2395.15    1266.05  -1.892   0.0718 .  
## P.1   229.49      13.71  16.741 5.28e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 827.2 on 22 degrees of freedom
## Multiple R-squared:  0.9982, Adjusted R-squared:  0.9981 
## F-statistic:  6215 on 2 and 22 DF,  p-value: < 2.2e-16

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

## List of 6
##  $ datos      :'data.frame': 24 obs. of  4 variables:
##   ..$ Y  : num [1:24] 12116 12655 13674 14202 15241 ...
##   ..$ X  : num [1:24] 66.5 68.2 70.8 73.9 77.2 ...
##   ..$ F  : num [1:24] 12181 12666 13436 14364 15269 ...
##   ..$ res: num [1:24] -65.3 -11.2 238.2 -161.5 -28.3 ...
##  $ Fregresores: num [1:24, 1:10] 1 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:24] "X1" "X2" "X3" "X4" ...
##   .. ..$ : chr [1:10] "C" "1" "2" "3" ...
##  $ Tregresores: num [1:24, 1:10] 0.204 0.204 0.204 0.204 0.204 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:10] "C" "1" "2" "3" ...
##  $ Nregresores: num 10
##  $ sse        : num 532799
##  $ gcv        : num 65241

plot(ts(E,frequency = 1, start = 1995),type="l",main="Consumo Energia Final (Ktep).España",ylab="")
lines(ts(fff.E$fitted, frequency = 1, start = 1995), type="l", col=2)
lines(ts(rdf.mod$datos$F, frequency = 1, start = 1995), type="l", col=3)
legend("top", ncol=3, c("CEF","Estimado FFF","Estimado RBS"), cex=0.6, bty="n", fill=c(1,2,3))

# REGRESION EN EL DOMINIO DE LA FRECUENCIA FILTRANDO LAS ALTAS FRECUENCIAS 
rdf2.mod=rdf2(E,P)

## Warning in rbind(C = c(1, rep(0, (n - 1))), cx, cx.2): number of columns of
## result is not a multiple of vector length (arg 1)

str(rdf2.mod)

## List of 5
##  $ datos      :'data.frame': 24 obs. of  4 variables:
##   ..$ Y  : num [1:24] 12116 12655 13674 14202 15241 ...
##   ..$ X  : num [1:24] 66.5 68.2 70.8 73.9 77.2 ...
##   ..$ F  : num [1:24] 12166 12698 13383 14455 15204 ...
##   ..$ res: num [1:24] -50.4 -43.4 291 -253.4 36.5 ...
##  $ Fregresores: num [1:24, 1:13] 1 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:13] "C" "" "" "" ...
##  $ Tregresores: num [1:24, 1:13] 0.2041 0.6273 0.0985 0.2887 -0.2959 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:13] "C" "" "" "" ...
##  $ Nregresores: num 13
##  $ Betas      : num [1:13, 1] 262.96 927.13 -40.19 1.21 -13.11 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:13] "C" "" "" "" ...
##   .. ..$ : NULL

plot(ts(E,frequency = 1, start = 1995),type="l",main="Consumo Energia Final (Ktep).España",ylab="")
lines(ts(rdf.mod$datos$F, frequency = 1, start = 1995), type="l", col=2)
lines(ts(rdf2.mod$datos$F, frequency = 1, start = 1995), type="l", col=3)
legend("top", ncol=3, c("CEF","Estimado RBS","Estimado RBS2"), cex=0.6, bty="n", fill=c(1,2,3))

Descomposición de una serie temporal en el dominio de la frecuencia.

El método de regresión de en el dominio de la frecuencia puede usarse para descomponer una serie de tiempo en componentes estacionales, de tendencia e irregulares de una serie de tiempo $y_t$ de frecuencia $b$ o número de veces en cada intervalo de tiempo unitario. Un valor de 7 para la frecuencia indicaría que los datos se muestrean diariamente y el período de tiempo natural es una semana, o 4 y 12 cuando los datos se muestrean trimestral y mensualmente y el período de tiempo natural es un año.

Si las observaciones se toman a intervalos iguales de longitud, $\bigtriangleup t$, entonces la frecuencia angular es $\omega = \frac {\pi} {\bigtriangleup t}$. La frecuencia equivalente expresada en ciclos por unidad de tiempo es $f =\frac {\omega} {2\pi} = \frac {1}{2} \bigtriangleup t$. Con solo una observación por año, $\omega =\pi$ radianes por año o $f =\frac {1}{2}$ ciclo por año (1 ciclo por dos años), la variación con una longitud de onda de un año tiene frecuencia $\omega=2\pi$ radianes por año o $f = 1$ ciclo por año.

Por ejemplo, en una serie de tiempo mensual de observación $N = 100 $, los ciclos estacionales o la longitud de onda de un año tienen una frecuencia $f = \frac {100} {12} = 8.33$ ciclos para 100 fechas. Si las series de tiempo son 8 años completos, la frecuencia menos estacional es 1 ciclo por año u 8 ciclos para 96 observaciones. Los números enteros multiplicados son $2 \frac {N} {12}$, $3 \frac {N} {12}$ …., y la longitud de onda mínima de un año tiene frecuencia son $f \leq \frac {N} {12}$ .

La serie de tiempo $y_t$ puede estimarse utilizando la RBS y una variable auxiliar o indice temporal $t$:

\[\dot y=W x_tI_nW^T\dot \beta + WI_nW^T\dot u\]

donde $x_t= a + bt$ siendo $t = (1,2,3, .... N) _N$.

Las frecuencias asociadas a la tendencia, o a los ciclos estacionales se pueden filtran para obtener la estimación RBS de las señales de tendencia y estacionales. Las primeros $2\frac {N}{12} - 1$ filas de la matriz en una serie de frecue4ncia mensual, se utilizan para estimar los coeficientes de Fourier correspondientes a ciclos de baja frecuencia, ciclos de tendencia y filas $2\frac {N}{12}$ y $2\frac {N} {12} + 1$ se utilizan para estimar los coeficientes de Fourier de 1 ciclo por año. El número entero se multiplica por las filas $6\frac {N}{12}$, $6\frac {N}{12} + 1$, $8 \frac {N}{12}$ se usar para estimar la serie estacional.

En la libreria R “descomponer” encontramos las funciones para realizar este tipo de descomposición temporal.

Función descomponer(y,frequency,type)

En base al modelo de regresión en el dominio de la frecuencia descompone una serie $y_t$ en los factores de tendencia $TD$, estacionales $ST$, e irregulares $IR$.

La función se desarrolla en los siguientes pasos:

Se calcula el periodograma de la serie, y se ordena según el vector de frecuencias para crear diferentes indices de orden.
Se obtiene un modelo de tendencia, a partir de las frecuencias mayores que $\frac{n}{2*frequency}$ si la serie es par y las mayores que $\frac{n-1}{2*frequency}$ si la serie es impar. Para ello, se realiza la regresión en domininio de la frecuencia entre la serie $y_t$ y los regresores que se obtienen con la matriz auxiliar $Wx_tI_nW^T$, donde $x_t$ es el resultado de ajustar un modelo lineal del tipo $y_t=a+bt+e_t$ a la serie de datos (tipo=1) ó un modelo cuadrático del tipo $y_t=a+bt+ct^2+e_t$, en donde solo se consideran los regresores correspondientes a las diferentes frecuencias seleccionadas. Una vez obtenidos los parámetros del modelo, se calcula la serie en el dominio de la frecuencia que una vez convierten al dominio del tiempo da como resultado la serie de tendencia $TD$.
Se obtiene la serie residual $IRST=y_t-TD$, se y sobre esa serie se realiza una nueva selección de frecuencias, las correspondientes a los factores estacionales es decir:$\frac{n}{2*frequency}$, $\frac{2n}{2*frequency}$,$\frac{3n}{2*frequency}$, etc….. Se realiza la regresión en el dominio de la frecuencia entre $IRST$ y los regresores correspondientes a las frecuencias seleccionadas obtenidas a partir de a matriz auxiliar $Wx_tI_nW^T$, donde $x_t$ es el resultado de ajustar un modelo lineal del tipo $IRST=a+bt+e_t$ a la serie de datos (tipo=1) ? un modelo cuadr?tico del tipo $IRST=a+bt+ct^2+e_t$. Una vez obtenidos los parámetros del modelo, se calcula la serie en el dominio de la frecuencia que una vez convierten al dominio del tiempo da como resultado la serie de tendencia $ST$.
Se obtiene la serie irregular a partir de $IR=IRST-ST$.

y=BJsales.lead
frequency=4
type=2
descomponer <- function (y,frequency,type) {
  # Author: Francisco Parra Rodriguez
  # http://rpubs.com/PacoParra/24432
  # date:"y", frequency:"frequency". 
  # Use 7 for frequency when the data are sampled daily, and the natural time period is a week, 
  # or 4 and 12 when the data are sampled quarterly and monthly and the natural time period is a year.
  n <- length(y) 
  y <- matrix(y,ncol=1)
  M <- MW(n) #crea la matriz de harvey para los n datos 
  f1 <- NULL
  if(n%%2==0) {f2 <- n/(2*frequency)} else {
    f2 <- (n-1)/(2*frequency)}
  #Modelo para obtener serie con tendencia
  c <- seq(from=2, to=(2+(n/frequency) ))
   i <- seq(1:n)
  i2 <- i*i  
   i <- seq(1:n)
  i2 <- i*i  
  if (type==1)
  {eq <- lm(y~i)  
   z <- eq$fitted} else {
     if (type==2) eq <- lm(y~i+i2) 
     z <- eq$fitted} 
     cx <- M%*%diag(z)
     cx <- cx%*%t(M)
    id <- seq(1,n)
  S1 <- data.frame(cx)
  S2 <- S1[1:(2+(n/frequency)),]
  X <- as.matrix(S2)
  cy <- M%*%y
  B <- solve(X%*%t(X))%*%(X%*%cy)
  Y <- t(X)%*%B
  BTD <- B
  XTD <- t(M)%*%t(X)
  TD <- t(M)%*%Y
  # Genero la serie residual
  IRST <- y-TD
  # Realizo la regresion dependiente de la frecuenca utilizando como explicativa IRST.
  # modelo para obtener serie con  estacionalidad con trunc.
  frecuencia <- seq(0:(n/2)) 
  frecuencia <- frecuencia-1
  S <- data.frame(f1=frecuencia)
  sel <- subset(S,f1==trunc(2*f2))
  c <- seq(from=2,to=(n/f2))
  for (i in c) {sel1 <- subset(S,f1==i*trunc(2*f2))
                sel <- rbind(sel,sel1)}
  m1 <- c(sel$f1 * 2)
  m2 <- c(m1+1)
  c <- c(m1,m2)
  n3 <- length(c)
  d <- rep(1,n3)
  s <- data.frame(c,d)
  S=s[with(s, order(c)), ]
  l <- frequency*trunc(n/frequency)
  ML <- MW(l)
  i <- seq(1:l)
  i2 <- i*i  
  if (type==1)
   {eq <- lm(y[1:l]~i)  
   z <- eq$fitted} else {
     if (type==2) eq <- lm(y[1:l]~i+i2) 
     z <- eq$fitted}   
   cx <- ML%*%diag(z) #problema
   cx <- cx%*%t(ML)
  id <- seq(1,l)
  S1 <- data.frame(cx,c=id)
  S2 <- merge(S,S1,by.x="c",by.y="c")
  S3 <- rbind(c(1,1,cx[1,]),S2) 
  m <- l+2
  X1 <- S3[,3:m]
  # matriz de regresores a l
  X1 <- as.matrix(X1)
  # la paso al dominio del tiempo
  X2 <- data.frame(t(ML)%*%t(X1))
  if (n==l) X3 <- X2 else
  X3 <- rbind(X2,X2[1:(n-l),])
  # la paso al dominio de la frecuencia
  X4 <-M%*%as.matrix(X3)
  cy <- M%*%IRST
  B1 <- solve(t(X4)%*%X4)%*%(t(X4)%*%cy)
  Y <- X4%*%B1
  BST <- B1
  XST <- M%*%X4
  ST <- t(M)%*%Y
  TDST <- TD+ST
  IR <- IRST-ST  
  data <- data.frame(y,TDST,TD,ST,IR)
  regresoresTD <- data.frame(XTD)
  regresoresST <- data.frame(XST)
 list(datos=data,regresoresTD=regresoresTD,regresoresST=regresoresST,coeficientesTD=BTD,coeficientesST=BST)
  }

Función gdescomponer (y,freq,type,year,q))

Gráfico de la función descomponer.

gdescomponer <- function(y,freq,type,year,q) {
  # Author: Francisco Parra Rodriguez
  # http://econometria.wordpress.com/2013/08/21/estimation-of-time-varying-regression-coefficients/ 
  serie <- descomponer (y,freq,type)
  TdsT <- c(serie$datos$TDST)
  Td <- c(serie$datos$TD)
  sT <- c(serie$datos$ST)
  TDST <- ts(TdsT,frequency=freq,start = c(year,q))
  TD <- ts(Td,frequency=freq,start = c(year,q))
  ST <- ts(sT,frequency=freq,start = c(year,q))
  par(mfrow=c(3,1))
  plot (TDST)
  plot (TD)
  plot (ST)
}

Ejemplo 14

Realizamos una descomposicion temporal para el IPI de Cantabria, con una asociacion quadratica.

library(descomponer)
data(ipi)
datos <- descomponer(ipi,12,2)
plot(ts(datos$datos,frequency=12))

gdescomponer(ipi,12,1,2002,1)

Bibliografia

Alvarez N. (1990). Una aproximación a los ciclos económicos. Economía aplicada cuantitativa Cuadernos de la UNED. 1990: 143-201.

Artis M.; Clavel J.G.; Hoffmann M. y Nachane D. (2007) “Harmonic Regression Models: A Comparative Review with Applications”. Institute for Empirical Research in Economics”. University of Zurich. Working Paper Series .September 2007 SSRN-Harmonic Regression Models: A Comparative Review with Applications by Michael Artis; José Clavel; Mathias Hoffmann; DILIP NACHANE

Ashley, Richard A. (1984), “A Simple Test for Regression Parameter Instability,” Economic Inquiry 22, No. 2, 253-267.

Barnett; W. A. (1983). “New Indices of Money Supply and the Flexible Laurent Demand System.” J. Bus. and Econ. Statist. 7-23.

Bartlett M.S. (1954); "Problémes de l’analyse spectrale des séries temporelles stationnaries“. Inst. Statist. Univ. Paris 3; 119-34

Bartlett M.S. (1966). “An introduction to Stochastic Processes”; 2 end.Ed. Cambridge University Press.

Berndt; E. R.; and M. S. Khaled (1979). “Parametric Productivity Measurement and Choice Among Flexible Functional Forms.” J. Polit. Econ. 87.1220-45.

Bhansali; R.J. (1979) “A Mixed Spectrum Analysis of the Lynx Data;”Journal of the Royal Statistical Society; Ser. A; 142(2); 199-209.

Brunk; H.D. (1962). “On the range of the difference between hypothetical distribution function and Pyke’s modified empirical distribution function”. Ann. Marh. Statist. 33; 525-32.

Campbell; M.J. and Walker; A.M. (1977); A Survey of Statistical Work on the Mackenzie River Series of Annual Canadian Lynx Trappings for the Years 1821-1934 and a New Analysis;" Journal of the Royal Statistical Society; Ser. A; 140(4); 411-431.

Caves; W.; Christensen; L. R. y Tretheway; M. W. (1980): “Flexible cost functions for multiproduct firms”. The Review of Economics and Statistics. Nº 62; págs 477-481.

Contreras; D y Escolano J (1984): EI análisis espectral como instrumento para detectar la estacionalidad. ESTADISTICA ESPAÑOLA Núm. 104; 1984; págs. 101 a 144

Chalfant J. y Gallant A. (1985); “Estimating substitution elasticities with the Fourier costs function; some Monte Carlo results”; Journal of Econometrics; 28; 205-222.

Christensen; L. R.; D. W. Jorgenson; and L. J. Lau. “Transcendental Logarithmic Production Frontiers.” Rev. Econ. and Statist. 55(1973):28-45.

De Janvry; A. (1972) “The Generalized Power Function.” Amer J. Agr.Econ. 54. 234-37.

Deaton A. y Muellbauer J. (1980); “An Almost Ideal Demand System”; American Economic Review; 70; 312-326

Denny; M. (1974). “The Relationship Between Forms of the Production Function.” Can. J. Econ. 7.21-31.

Despotakis; K. A. (1986). “Economic Performance of Flexible Functional Forms.” Eur. Econ. Rev. 30.1107-43.

Diewert; W. E. (1971). “An Application of the Shephard Duality Theorem: A Generalized Leontief Production Function.” J. Polit. Econ. 79. 481-507.

Diewert W. y Wales T. (1987); “Flexible Functional Forms and Global Curvature Conditions”; conometrica; 55; No. 1; 43-68.

Diewert W. y Wales T. (1988); “A Normalized Quadratic Semiflexible Functional Form”; Journal of Econometrics; 37; 327-42.

Durbin J. (1969); “Test for serial correlation in regresion analisys based on the periodogram of least squared residuals”. Biometrica; vol 56; nº1; (mar-1969); pp 1-15.

Engle R. (1972). “Band spectrum regression”. Working Paper Department of economics. MIT. Nº26. November 1972.

Engle, Robert F. (1974), Band Spectrum Regression,International Economic Review 15,1-11.

Elbadawi; I.; A. R. Gallant; and G. Souza. (1983). “An Elasticity Can Be Estimated Consistently without A Priori Knowledge of Functional Form.” Econometrica 51.1731-51.

Fleissig A.; Kastns T. y Terrell D. (1997); “Semi-nonparametric estimates of substitution elasticities”; Economics Letters; 54; No. 3; 209-219.

Fuss; M.; D. McFadden; and Y. Mundlak. (1978). “A Survey of Functional Forms in the Economic Analysis of Production.” Production Economics: A Dual Approach toTheory and Applications; ed. M. Fuss and D. Mc-Fadden; pp. 219-68. Amsterdam: North-Holland Publishing Co.

Gallant; A. R.(1981) “On the Bias in Flexible Functional Forms and an Essentially Unbiased Form.” J. Econometrics 15(1981):211-45.

Gallant; A. R. (1982). “Unbiased Determination of Production Technologies.” J. Econometrics 20. 285-323.

Gallant; A. R.(1984) “The Fourier Flexible Form.” Amer. J. Agr. Econ.66(1984):204-15.

Gallant A. y Golub G. (1984); “Imposing Curvature Restrictions on Flexible Functional Forms”; Journal of Econometrics; 26; 295-321

Gonzalez G: Series de Fourier; Transformadas de Fourier y Aplicaciones; disponible en www.dmae.upm.es/WebpersonalBartolo/VariableCompleja/VCParteI/9_Series_de_Fourier.ppt

Griffin R.C.; Montgomery J.M.; and Rister M.E. (1987). Selecting Functional Form in Production Function. Análisis Western Journal of Agricultural Economics; 12(2): 216-227

Guilkey D. y Lovell C. (1980); “On the Flexibility of the Translog Approximation”; International Economic Review; 21; 137-147.

Halter; A. N.; H. O. Carter; and J. G. Hocking. (1957) “A Note on the Transcendental Production Function.” J. Farm Econ. 39. 966-74.

Hannan; E. J. (1963); “Regression for Time Series;” in Time Series Analysis ed. by M. Rosenblatt; John Wiley; 1963.

Hannan; E.J.(1971); “Non-linear Time Series Regression;” Journal ofApplied Probability; 8; 767-780.

Hannan; E.J.(1973); “The Asymptotic Theory of Linear Time Series Models;” Journal of Applied Probability; 10; 130-145.

Harvey; A.C. (1978): “Linear Regression in Frequency Domain”. International Economic Review 19(2), 507-512.

Heady; E. 0.; and J. L. Dillon; eds. (1961). “Agricultural Production Functions”. Ames: Iowa State University Press.

Hui;T and Ashley R. (1999): “An elementary method for detecting and modelling regression parameter variation across frequencies whit an application to testing the permanent income hypothesis”. Macroeconomic Dynamics / Volume 3 / Issue 01 / March 1999, pp 69-831999 Cambridge University Press.(http://ashleymac.econ.vt.edu/working_papers/e9702.pdf)

Kay; S.M. and Marple; S.L. (1981); “Spectrum Analysis - A Modern Perspective;” Proceedings of the IEEE; 69; 1380-1419.

Ladd; G. W. (1979). “Artistic Tools for Scientific Minds.” Amer.J. Agr. Econ. 61.1-11.

Lau; L. J. (1978) “Testing and Imposing Monotonicity; Convexity and Quasi-Convexity Constraints.” Production Economics: A Dual Approach to Theory and Applications; ed. M. Fuss and D. McFadden; pp. 409-53. Amsterdam: North-Holland Publishing Co.

Lilyan E. Fulginiti; Richard K. Perrin and Bingxin Yu (2003): “Institution and agricultural productivity in sub-sahara Africa”. 25 International Conference of Agricultural Economists (IAAE). Durban South Africa.

Marple; S. (1987) “Digital Spectral Analysis with Applications”; Prentice-Hall; Englewood Clios; N.J.

Martínez M., Gómez L., Serrano J.A., Vila J., Gómez J. (2009): “Filtros Digitales”. Curso 2009-2010. Escola Tècnica Superior Ingenieria. Departament d´Enginyeria Electrònica Universidad de Valencia.

Melis F. (1991): “La Estimación del ritmo de variación de las series económicas”. Estadística Española. Vol 22. Num. 126; 1991.

McFadden; D. (1963). “Constant Elasticity of Substitution Production Functions.” Rev. Econ. Stud. 30.73-83. Maddala; G. S. Econometrics. New York: McGraw-Hill Book Co.; 1977.

Parra F (2014): Seasonal Adjustment by Frequency Analysis. Package R Version 1.1. URL:http://cran.r-project.org/web/packages/descomponer/index.html

Pyke R. (1959). “The supremun and infimum of the Poisson process”. Ann. Math. Statist.30; 568-76.

Pisarenko; V.F. (1973); “The Retrieval of Harmonics from a Covariance Function;” Geophysical Journal of Royal Astronomical Society; 33; 347-366.

Priestley; M.B. (1964); “Estimation of the Spectral Density Function in the Presence of Harmonic Components;” Journal of the Royal Statistical Society; Ser. B; 26; 123-132.

Priestley; M.B. (1981); “Spectral Analysis and Time Series; Academic Press; London.

Ramajo J. “Avances recientes en el análisis econométrico de la demanda “.http://eco.unex.es/jramajo/ivcnea.pdf

Silberberg; E. (1978). The Structure of Economics: A Mathematical Analysis. New York: McGraw-Hill Book Co.

Truong-Van; B. (1990); “A New Approach to Frequency Analysis with Amplied armonics;” Journal of the Royal Statistical Society; Ser.B; 52(1); 203-221.

Tweeten; L. (1983). “Hypothesis Testing in Economic Science.” Amer. J.Agr. Econ. 65.548-52.

Uzawa; H. (1962). “Production Functions with Constant Elasticity of Substitution.” Rev. Econ. Stud. 29.291-99.

Venables and Ripley (2002), “Modern Applied Statistics with S” (4th edition, 2002).

Villarreal F. (2005): “Elementos teóricos del ajuste estacional de series económicas utilizando X-12-ARIMA y TRAMO-SEATS”. Serie 38. Estudios estadísticos y prospectivos División de Estadística y Proyecciones Económicas. Santiago de Chile; Diciembre del 2005. www.eclac.org/publicaciones/xml/9/24099/lcl2457e.pdf

Whittle; P. (1952); “Tests of Fit in Time Series;” Biometrika; 39; 309-318.

White; H. (1980). “Using Least Squares to Approximate Unknown Regression Functions.” Int. Econ. Rev. 21.149-70.

Young; P.C.; Pedregal; D.J. and W.Tych (1999); “Dynamic harmonic regression;” Journal of Forecasting; 18; 369-394

Econometria con Series de Fourier

Francisco Parra

14/10/2020