Descomposicion temporal

by Francisco Parra

25 Marzo 2015

1. Datos

Series:

IPI. Base 1005. Enero 2002-Septiembre 2013

tasa interanual de crecimiento del PIB de cantabria en indices de volumen (pibiv). 1T 2001 - 4T 2013

Población (pob). EPA 1T 2005 - 1T 2014

Activos (acti). EPA 1T 2005 - 1T 2014

Ocupados (ocu). EPA 1T 2005 - 1T 2014

Parados (par). EPA 1T 2005 - 1T 2014

Parados primer empleo (par1). EPA 1T 2005 - 1T 2014

Inactivos (ina). EPA 1T 2005 - 1T 2014

Afiliados a todos los regímenes (trg). Enero 1999-Mayo 2014

Afiliados régimen general (tg). Enero 1999-Mayo 2014

Afiliados régimen autónomo (reta). Enero 1999-Mayo 2014

IBEX-35 Datos diarios.

ipi <- read.csv("ipi.csv", header = FALSE, 
    sep = ",", dec = ".")
ipi <- ipi$V1
pibiv <- read.csv("pibiv.csv", header = FALSE, 
    sep = " ", dec = ",")
pib <- pibiv$V1
epa <- read.csv("epa.csv", header = FALSE, 
    sep = ";", dec = ",")
pob <- epa$V1
act <- epa$V2
ocu <- epa$V3
par <- epa$V4
par1 <- epa$V5
ina <- epa$V6
afi <- read.csv("afi.csv", header = FALSE, 
    sep = ";", dec = ".")
trg<- afi$V1
rg <- afi$V2
reta <- afi$V3
a <- rnorm(2000,0,2)
b <- seq(1:2000)
c <- a+ 0.1*b-0.00002*b^2
ibex <- read.csv("ibex.csv", header=FALSE, sep = ",", dec = ".")
ibex <- ibex$V1[1:2000]

2. Representación datos

## plot time series
plot(ts(ipi,frequency = 12), type="l")

plot(ts(pib,frequency = 4), type="l")

plot(ts(pob,frequency = 4), type="l")

plot(ts(act,frequency = 4), type="l")

plot(ts(ocu,frequency = 4), type="l")

plot(ts(par,frequency = 4), type="l")

plot(ts(par1,frequency = 4), type="l")

plot(ts(ina,frequency = 4), type="l")

plot(ts(trg,frequency = 12), type="l")

plot(ts(rg,frequency = 12), type="l")

plot(ts(reta,frequency = 12), type="l")

plot(ts(ibex,frequency=365),  type="l")

2. Funciones R

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)

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

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

donde \(\eta\) es la media de la serie, \(a_j\) y \(b_j\) son su amplitud,\(f\) son las frecuencias que del conjunto de las \(n\) observaciones,\(t\) es un indice de tiempo que va de 1 a N, siendo N el numero de periodos para los cuales tenemos observaciones en el conjunto de datos, el cociente \(\frac{ft}n)\) convierte cada valor de \(t\) en escala de tiempo en proporciones de \(2n\) y rango \(j\) desde \(1\) hasta \(n\) siendo \(n=\frac{N}2\) (es decir, 0,5 ciclos por intervalo de tiempo). Las dinámica de las altas frecuencias (los valores más altos de f) corresponden a los ciclos cortos en tanto que la dinámica de la bajas frecuencias (pequeños valores de f) van a corresponder con los ciclos largos. Si nosotros hacemos que \(\frac{ft}n=w\) la ecuación (1) quedaría, asi :

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

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

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

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

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

Matriz \(W\)

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

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

Función cdf(a)

Otiene 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) La multiplicación de dos series armónicas de diferente frecuencia:

\([a_j\cos (\omega_j)+b_j\sin (\omega_j)]x [a_i\cos (\omega_i)+b_i\sin (\omega_i)]\)

da como resultado la siguiente suma: \[\begin{equation} \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} \end{equation}\]

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

\[\begin{equation} \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} \end{equation}\]

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{equation} \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} \end{equation}\]

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

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

Función periodograma (a)

Calcula y presenta el espectro de la serie “a”

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

\(\hat a= Wa\)

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

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

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

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

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 gperiodogrma (a)

Presenta gráficamente el espectro de la variabe 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")}

Función que desconpone una serie temporal en tendencia, estacionalidad y componente irregular

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 ó mayeros 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\).

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)
  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) ))
  #Use the "sort.data.frame" function, Kevin Wright. Package taRifx
  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 <- cdf(z)
  id <- seq(1,n)
  S1 <- data.frame(cx)
  S2 <- S1[1:(2+(n/frequency)),]
  X <- as.matrix(S2)
  cy <- gdf(y)
  B <- solve(X%*%t(X))%*%(X%*%cy)
  Y <- t(X)%*%B
  BTD <- B
  XTD <- t(MW(n))%*%t(X)
  TD <- gdt(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 <- sort.data.frame (s,formula=~c)
  #Use the "sort.data.frame" function, Kevin Wright. Package taRifx
  # Se realiza el ejercicio con los datos del año completo
  l <- frequency*trunc(n/frequency)
  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 <- cdf(z)
  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(MW(l))%*%t(X1))
  if (n==l) X3 <- X2 else
    X3 <- rbind(X2,X2[1:(n-l),])
  # la paso al dominio de la frecuencia
  X4 <-MW(n)%*%as.matrix(X3)
  cy <- gdf(IRST)
  B1 <- solve(t(X4)%*%X4)%*%(t(X4)%*%cy)
  Y <- X4%*%B1
  BST <- B1
  XST <- t(MW(n))%*%X4
  ST <- gdt(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)}

4. Resultados y comparación con “descomponse” y “stl”

library(descomponer)

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

## Loading required package: taRifx

## Warning: package 'taRifx' was built under R version 3.2.0

## 
## Attaching package: 'descomponer'
## 
## The following objects are masked _by_ '.GlobalEnv':
## 
##     cdf, descomponer, gdf, gdt, gperiodograma, MW, periodograma

par(mfrow=c(5,1))
plot(ts(descomponer(ipi,12,1)$datos,frequency = 12))

ipi <- ts(ipi,frequency = 12)
plot(decompose(ipi,type="additive"))

plot(stl(ipi, s.window=12))

plot(ts(descomponer(pib,4,1)$datos,frequency = 4))

pib <- ts(pib,frequency = 4)
plot(decompose(pib,type="additive"))

plot(stl(pib, s.window=4))

plot(ts(descomponer(pob,4,1)$datos,frequency = 4))

pob <- ts(pob,frequency = 4)
plot(decompose(pob,type="additive"))

plot(stl(pob, s.window=4))

plot(ts(descomponer(act,4,1)$datos,frequency = 4))

act <- ts(act,frequency = 4)
plot(decompose(act,type="additive"))

plot(stl(act, s.window=4))

plot(ts(descomponer(ocu,4,1)$datos,frequency = 4))

ocu <- ts(ocu,frequency = 4)
plot(decompose(ocu,type="additive"))

plot(stl(ocu, s.window=4))

plot(ts(descomponer(par,4,1)$datos,frequency = 4))

par <- ts(par,frequency = 4)
plot(decompose(par,type="additive"))

plot(stl(par, s.window=4))

plot(ts(descomponer(par1,4,1)$datos,frequency = 4))

par1 <- ts(par1,frequency = 4)
plot(decompose(par1,type="additive"))

plot(stl(par1, s.window=4))

plot(ts(descomponer(ina,4,1)$datos,frequency = 4))

ina <- ts(ina,frequency = 4)
plot(decompose(ina,type="additive"))

plot(stl(ina, s.window=4))

plot(ts(descomponer(trg,12,1)$datos,frequency = 12))

trg <- ts(trg,frequency = 12)
plot(decompose(trg,type="additive"))

plot(stl(trg, s.window=12))

plot(ts(descomponer(rg,12,1)$datos,frequency = 12))

rg <- ts(rg,frequency = 12)
plot(decompose(rg,type="additive"))

plot(stl(rg, s.window=12))

plot(ts(descomponer(reta,12,1)$datos,frequency = 12))

reta <- ts(reta,frequency = 12)
plot(decompose(reta,type="additive"))

plot(stl(reta, s.window=12))

plot(ts(descomponer(ibex,365,1)$datos,frequency=365))

ibex <- ts(ibex,frequency = 365)
plot(decompose(ibex,type="additive"))

plot(stl(ibex, s.window=365))

Bibliografía

Cleveland R.B,Cleveland W.S, McRae J.E. (1990)“STL: A Seasonal-Trend Decomposition Procedure Based on Loess”. Journal of Oficial Statistics. Vol 6 nº1 1990 pp-3-75. ( http://cs.wellesley.edu/~cs315/Papers/stl%20statistical%20model.pdf )

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

Hannan, E.J. (1963), Regression for Time Series, in Rosenblatt, M. (ed.), Time Series Analysis, New York, John Wiley.

Harvey, A.C. (1978), Linear Regression in the Frequency Domain, International Economic Review, 19, 507-512.

Parra, F. (2013), Regresión con paramétros dependientes del tiempo, (http://econometria.wordpress.com/2013/07/29/estimacion-con-parametros-dependientes-del-tiempo/)

```