Tema 4: Regresión con series de tiempo

Curso: Análisis de series temporales

Paquetes de R

Para este tema, se necesita cargar estos paquetes:

library(ggplot2)
library(forecast)
library(fpp2)
library(astsa)
library(car)
library(astsa)
library(dplyr)
library(gamair)
library(mgcv)

Introducción

Contenido

  1. Introducción

  2. Modelos de regresión

  3. Regresión espuria

  4. El tiempo como covariables

  5. Pronósticos usando regresión

  6. Extensiones de modelos de regresión

Introducción

  • La idea es ajustar un modelo de regresión para la serie temporal \(Y_t, t=1,...,T\) utlizando un conjunto de \(p\) covariables: \(X_1,...,X_p\).

  • Por ejemplo:

    • \(Y\) una serie mensual de ventas con la variable independiente \(X_1\) gasto mensual en anuncios.
    • \(Y\) una serie diaria de demanda de energía eléctrica con las variables independientes: \(X_1\) temperatura y \(X_2\) el día de la semana.

Modelos de regresión

Contenido

  1. Introducción

  2. Modelos de regresión

  3. Regresión espuria

  4. El tiempo como covariables

  5. Pronósticos usando regresión

  6. Extensiones de modelos de regresión

Regresión lineal simple

Un modelo de regresión lineal simple establece una relación lineal entre una variable dependiente \(Y\) y una sola variable predictora \(X\):

\[ Y_t = \beta_0 + \beta_1 X_t + \epsilon_t, \] donde los coeficientes \(\beta_0\) y \(\beta_1\) denotan la intersección y la pendientes, respectivamente;
\(\varepsilon_t \overset{\text{iid}}{\sim} N(0,\sigma^2)\).

  • La intercepción \(\beta_0\) representa el valor predicho de \(y\) cuando \(X=0\).
  • La pendiente \(\beta_1\) representa el cambio promedio previsto en \(Y\) resultante de un aumento de una unidad en \(X\).

Ejemplo simulado de un modelo de regresión lineal simple (Fig 5.1 de Hyndman)

Ejemplo

  • La serie de cambios porcentuales trimestrales (tasas de crecimiento) del gasto de consumo personal real (\(Y_t\), variable dependiente), y el cambio porcentual de los ingresos disponibles(\(X_t\), covariable), para EE.UU. desde 1970 a 2016.
head(uschange)
        Consumption     Income Production   Savings Unemployment
1970 Q1   0.6159862  0.9722610 -2.4527003 4.8103115          0.9
1970 Q2   0.4603757  1.1690847 -0.5515251 7.2879923          0.5
1970 Q3   0.8767914  1.5532705 -0.3587079 7.2890131          0.5
1970 Q4  -0.2742451 -0.2552724 -2.1854549 0.9852296          0.7
1971 Q1   1.8973708  1.9871536  1.9097341 3.6577706         -0.1
1971 Q2   0.9119929  1.4473342  0.9015358 6.0513418         -0.1
autoplot(uschange[,c("Consumption","Income")]) +
  ylab("% change") + xlab("Year")

uschange %>%
  as.data.frame() %>%
  ggplot(aes(x=Income, y=Consumption)) +
    ylab("Consumption (quarterly % change)") +
    xlab("Income (quarterly % change)") +
    geom_point() +
    geom_smooth(method="lm", se=FALSE)

  • Más adelante vemos que tslm es programado para series temporales, pero el cálculo es lo mismo.
lm(Consumption ~ Income, data=uschange)

Call:
lm(formula = Consumption ~ Income, data = uschange)

Coefficients:
(Intercept)       Income  
     0.5451       0.2806  
tslm(Consumption ~ Income, data=uschange)

Call:
tslm(formula = Consumption ~ Income, data = uschange)

Coefficients:
(Intercept)       Income  
     0.5451       0.2806  

\[Y_t=0.55 + 0.28X_t + \epsilon_t\]

  • El coeficiente de pendiente muestra que un aumento de una unidad en \(X\)(un aumento de 1 punto porcentual en el ingreso personal disponible) resulta en un promedio de 0.28 unidades de aumento en \(Y\).

Regresión lineal múltiple

  • La forma general de un modelo de regresión lineal múltiple es:

\[Y_t=\beta_0+\beta_1 X_{t,1}+\beta_2 X_{t,2}+...+\beta_p X_{t,p}+\epsilon_t, t=1,...,T,\]

donde \(Y_t\) es la variable a pronosticar y \(X_1,...,X_p\) son los \(p\) variables predictoras. Las variables predictoras pueden ser numéricas o categóricas (con el manejo apropiado de factores).
Los coeficientes \(\beta_1,...,\beta_p\) miden el efecto de cada predictor después de tener en cuenta los efectos de todos los demás predictores del modelo. Por lo tanto, los coeficientes miden los efectos marginales de las variables predictoras.

  • El modelo de regresión lineal múltiple en su forma matricial:

\[ Y=X \beta+\epsilon \] donde

\[Y=\left[ \begin{array}{c}Y_1 \\ \vdots \\Y_T \end{array} \right],~~ X= \left(\begin{array}{ccccc} 1& X_{11}& X_{12} & ... & X_{1p}\\ 1 & X_{21}& X_{22} & ... &X_{2p}\\ \vdots& \vdots & \ddots &\vdots& \vdots\\ 1& X_{T1}& X_{T2} & ... &X_{Tp} \end{array}\right),\] \[\beta=\left[ \begin{array}{c}\beta_0 \\ \vdots \\\beta_T \end{array} \right],~~\epsilon=\left[ \begin{array}{c}\epsilon_1 \\ \vdots \\\epsilon_T \end{array} \right].\]

Supocisiones del modelo:

  • La relación entre la variable de pronóstico y las variables predictoras satisface esta ecuación lineal.

  • Los errores \(\varepsilon_1,...,\varepsilon_T\):

    • tienen media cero,
    • no están autocorrelacionados,
    • no están relacionados con las variables predictoras
  • Los errores se distribuyan normalmente con una varianza constante \(\sigma^2\).

  • Cada predictor \(X_i, i=1,...,p\) supone que es observado y fijo, i.e. no es una variable aleatoria.

Tópicos importantes:

  • Estimación:
    • Por mínimos cuadrados.
    • Por máxima verosimilitud.
  • Selección de variables
  • Diagnósticos
  • Medidas remediales
  • Estimación por mínimos cuadrados: minimizar \[\sum_{t=1}^T \epsilon_t^2=\sum_{t=1}^T [y_t-(\beta_0+\beta_1 x_{1,t}+...+\beta_k x_{k,t})]^2,\] en función de \(\beta_0,...\beta_k\).

  • Como resultado: \[\hat{\beta}=(X^\top X)^{-1}X^\top Y.\]

  • El estimador de máxima verosimilitud es equivalente.

Ejemplo

  • El objetivo es generar pronósticos del consumo más precisos usando otros predictores, además del ingreso personal.
autoplot(uschange[,3:5], facets = TRUE, colour=TRUE) +
  ylab("") + xlab("Year") +
  guides(colour="none")

Variaciones porcentuales trimestrales en la producción industrial y ahorros personales y variaciones trimestrales en la tasa de desempleo de los EE. UU. Durante el período 1970T1-2016T3 (Hyndman)

uschange %>%
  as.data.frame() %>%
  GGally::ggpairs()

Matriz de diagrama de dispersión del gasto de consumo de EE. UU. y los cuatro predictores.

modreg <- tslm(
  Consumption ~ Income + Production + Unemployment + Savings,
  data=uschange)
summary(modreg)

Call:
tslm(formula = Consumption ~ Income + Production + Unemployment + 
    Savings, data = uschange)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.88296 -0.17638 -0.03679  0.15251  1.20553 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.26729    0.03721   7.184 1.68e-11 ***
Income        0.71449    0.04219  16.934  < 2e-16 ***
Production    0.04589    0.02588   1.773   0.0778 .  
Unemployment -0.20477    0.10550  -1.941   0.0538 .  
Savings      -0.04527    0.00278 -16.287  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3286 on 182 degrees of freedom
Multiple R-squared:  0.754, Adjusted R-squared:  0.7486 
F-statistic: 139.5 on 4 and 182 DF,  p-value: < 2.2e-16
autoplot(uschange[,'Consumption'], series="Data") +
  forecast::autolayer(fitted(modreg), series="Fitted") +
  xlab("Year") + ylab("") +
  ggtitle("Cambio porcentual de gastos de consumo en EE.UUU ") +
  guides(colour=guide_legend(title=" "))
checkresiduals(modreg)

    Breusch-Godfrey test for serial correlation of order up to 8

data:  Residuals from Linear regression model
LM test = 14.874, df = 8, p-value = 0.06163

Regresión espuria

Contenido

  1. Introducción

  2. Modelos de regresión

  3. Regresión espuria

  4. El tiempo como covariables

  5. Pronósticos usando regresión

  6. Extensiones de modelos de regresión

Regresión espuria

  • Cuando las series no son estacionarias1, los resultados de regresión no son confiables y presenta lo que se llama correlación espuria.
  • Los datos de las series cronológicas de tendencias pueden parecer relacionados.

  • Por ejemplo, los pasajeros aéreos en Australia tienen una correlación positiva con la producción de arroz en Guinea.

fit <- tslm(aussies ~ guinearice)
summary(fit)

Call:
tslm(formula = aussies ~ guinearice)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.9448 -1.8917 -0.3272  1.8620 10.4210 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -7.493      1.203  -6.229 2.25e-07 ***
guinearice    40.288      1.337  30.135  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.239 on 40 degrees of freedom
Multiple R-squared:  0.9578,    Adjusted R-squared:  0.9568 
F-statistic: 908.1 on 1 and 40 DF,  p-value: < 2.2e-16
checkresiduals(fit)


    Breusch-Godfrey test for serial correlation of order up to 8

data:  Residuals from Linear regression model
LM test = 28.813, df = 8, p-value = 0.000342

El tiempo como covariables

Contenido

  1. Introducción

  2. Modelos de regresión

  3. Regresión espuria

  4. El tiempo como covariables

  5. Pronósticos usando regresión

  6. Extensiones de modelos de regresión

El tiempo como covariables

El tiempo también sirve como variable predictora:

  • Modelos de tendencia (lineales o no lineales)
  • Modelos estacionalidad (variables indicadoras)

Veremos más adelante en regresión dinámica: - Análisis de intervención. - Estrategias para valores extremos.

Modelos de tendencia

  • Como las variables independiente son asumidas como fijas, se puede utilizar el tiempo como una variable independiente.

  • Los modelos más frecuentes:

    • Tendencia lineal: \[Y_t=\beta_0+\beta_1 t + \epsilon_t\]
    • Tendencia cuadrática: \[Y_t=\beta_0+\beta_1 t +\beta_2 t^2 + \epsilon_t\]
  • Regresión no lineal.

    • Por ejemplo: LOESS.

Vamos a ajustar un modelo de tendencia cuadrática a la serie de graduados de ITCR de 1975 a 2002: \[Y_t=\beta_0+\beta_1 t +\beta_2 t^2 + \epsilon_t, t=1,...,T\]

itcrgrad<-read.csv("ITCR.csv",sep=",")
y<-ts(itcrgrad$graduados,start=1975)
tiempo<-seq(1,length(y))
tiempo2<-tiempo^2
mod<-lm(y~tiempo+tiempo2)

autoplot(y) +
  ylab("ITCR") +
  autolayer(ts(mod$fitted.values,start=1975), series = "ajustado") 

summary(mod)

Call:
lm(formula = y ~ tiempo + tiempo2)

Residuals:
     Min       1Q   Median       3Q      Max 
-117.932  -49.683    4.506   43.234  155.390 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 127.0021    43.6818   2.907  0.00753 ** 
tiempo      -12.7887     6.9427  -1.842  0.07736 .  
tiempo2       1.5612     0.2323   6.720 4.83e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 71.61 on 25 degrees of freedom
Multiple R-squared:  0.944, Adjusted R-squared:  0.9395 
F-statistic: 210.6 on 2 and 25 DF,  p-value: 2.268e-16
  • Supuestos:
    • Normalidad
    • Homoscedasticidad
    • Autocorrelación en el tiempo
hist(mod$residuals,main="",xlab="residuales")

car::qqPlot(mod$residuals,ylab="residuales")

[1]  8 21
shapiro.test(mod$residuals)

    Shapiro-Wilk normality test

data:  mod$residuals
W = 0.97509, p-value = 0.7211
  • Los residuales en el tiempo:
ts.plot(mod$residuals,ylab="residuales") 
checkresiduals(mod)

    Breusch-Godfrey test for serial correlation of order up to 6

data:  Residuals
LM test = 3.9215, df = 6, p-value = 0.6873

Transformación en modelos estacionales

  • Si la magnitud del cambio estacional se mantiene aproximadamente constante con el cambio del nivel de la serie, se dice que la variación estacional es constante (aditiva).
  • Si la variación estacional aumenta proporcionalmente con el nivel de la serie, se dice que la variación estacional es multiplicativa.

Transformaciones usadas en la práctica:

  • \(W_t=Y_t^\alpha,~~\text{con} -1<\alpha<1\)

  • \(W_t=\ln Y_t\)

Ejemplo de turistas

turistas<-read.csv("turistas.csv",sep=";")
y<-ts(turistas$turistas,start=c(1991,1),frequency=12)

Sin transformación \((Y_t)\)

autoplot(y) 

Con transformación \((\ln Y_t)\)

autoplot(log(y))  

Modelos de series estacionales mediante variables indicadoras

  • La idea es utilizar un factor (variables indicadoras) como predictor para indicar el periodo de la estacionalidad.
  • Ejemplo de turistas (datos mensuales y con tendencia cuadrática)

\[ln Y_t=\alpha_0+\alpha_1 t + \alpha_2 t^2 +\beta_1 I_{1}+...+\beta_11 I_{11}+\epsilon_t\]

w<-log(y)
tiempo<-seq(1,length(y))
tiempo2<-tiempo^2
mes<-rep(seq(1,12),10)
mes<-as.factor(mes)

datos1<-data.frame(w,tiempo,tiempo2,mes)
head(datos1)
         w tiempo tiempo2 mes
1 10.76223      1       1   1
2 10.72731      2       4   2
3 10.78324      3       9   3
4 10.50832      4      16   4
5 10.36113      5      25   5
6 10.51486      6      36   6
tail(datos1)
           w tiempo tiempo2 mes
115 11.42852    115   13225   7
116 11.26863    116   13456   8
117 11.08610    117   13689   9
118 11.13942    118   13924  10
119 11.45100    119   14161  11
120 11.68831    120   14400  12
mod<-lm(w~tiempo+tiempo2+mes,datos1)
round(summary(mod)$coefficients,4)
            Estimate Std. Error  t value Pr(>|t|)
(Intercept)  10.9420     0.0294 372.1819   0.0000
tiempo        0.0086     0.0008  11.1807   0.0000
tiempo2       0.0000     0.0000  -3.2291   0.0017
mes2         -0.1013     0.0324  -3.1246   0.0023
mes3         -0.1181     0.0324  -3.6408   0.0004
mes4         -0.3467     0.0324 -10.6900   0.0000
mes5         -0.5050     0.0324 -15.5680   0.0000
mes6         -0.4243     0.0324 -13.0771   0.0000
mes7         -0.2037     0.0325  -6.2784   0.0000
mes8         -0.3326     0.0325 -10.2482   0.0000
mes9         -0.6091     0.0325 -18.7601   0.0000
mes10        -0.5293     0.0325 -16.2990   0.0000
mes11        -0.3263     0.0325 -10.0431   0.0000
mes12        -0.1042     0.0325  -3.2068   0.0018
mod2<-tslm(w~trend+I(trend^2)+season,datos1)
round(summary(mod2)$coefficients,4)
            Estimate Std. Error  t value Pr(>|t|)
(Intercept)  10.9420     0.0294 372.1819   0.0000
trend         0.0086     0.0008  11.1807   0.0000
I(trend^2)    0.0000     0.0000  -3.2291   0.0017
season2      -0.1013     0.0324  -3.1246   0.0023
season3      -0.1181     0.0324  -3.6408   0.0004
season4      -0.3467     0.0324 -10.6900   0.0000
season5      -0.5050     0.0324 -15.5680   0.0000
season6      -0.4243     0.0324 -13.0771   0.0000
season7      -0.2037     0.0325  -6.2784   0.0000
season8      -0.3326     0.0325 -10.2482   0.0000
season9      -0.6091     0.0325 -18.7601   0.0000
season10     -0.5293     0.0325 -16.2990   0.0000
season11     -0.3263     0.0325 -10.0431   0.0000
season12     -0.1042     0.0325  -3.2068   0.0018
pronostico<-forecast(mod2,h=12)

autoplot(w) +
  ylab("turistas") +
  autolayer(mod2$fitted.values, series = "ajustado") +
  autolayer(pronostico, series = "pronostico")

Pronósticos usando regresión

Contenido

  1. Introducción

  2. Modelos de regresión

  3. Regresión espuria

  4. El tiempo como covariables

  5. Pronósticos usando regresión

  6. Extensiones de modelos de regresión

Pronósticos usando regresión

y<-ts(turistas$turistas,start=c(1991,1),frequency=12)
y.train<-window(y,start=c(1991,1),end=c(1999,12))
y.test<-window(y,start=c(2000,1),end=c(2000,12))
mod4<-tslm(y.train~trend+I(trend^2)+season)
summary(mod4)

Call:
tslm(formula = y.train ~ trend + I(trend^2) + season)

Residuals:
     Min       1Q   Median       3Q      Max 
-16096.7  -3411.3    833.9   3542.2  15635.4 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.285e+04  2.312e+03  27.188  < 2e-16 ***
trend        4.557e+02  6.746e+01   6.755 1.18e-09 ***
I(trend^2)  -5.924e-01  5.994e-01  -0.988 0.325522    
season2     -8.809e+03  2.552e+03  -3.452 0.000836 ***
season3     -1.006e+04  2.552e+03  -3.942 0.000156 ***
season4     -2.517e+04  2.552e+03  -9.860 3.65e-16 ***
season5     -3.402e+04  2.553e+03 -13.325  < 2e-16 ***
season6     -2.947e+04  2.553e+03 -11.540  < 2e-16 ***
season7     -1.574e+04  2.554e+03  -6.162 1.77e-08 ***
season8     -2.419e+04  2.555e+03  -9.470 2.46e-15 ***
season9     -3.947e+04  2.556e+03 -15.446  < 2e-16 ***
season10    -3.551e+04  2.556e+03 -13.890  < 2e-16 ***
season11    -2.476e+04  2.558e+03  -9.680 8.82e-16 ***
season12    -8.906e+03  2.559e+03  -3.481 0.000761 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5414 on 94 degrees of freedom
Multiple R-squared:  0.9158,    Adjusted R-squared:  0.9041 
F-statistic:  78.6 on 13 and 94 DF,  p-value: < 2.2e-16
pronostico<-forecast(mod4,h=12)
pronostico
         Point Forecast    Lo 80     Hi 80    Lo 95     Hi 95
Jan 2000      105476.42 97842.25 113110.58 93731.85 117220.99
Feb 2000       96993.44 89340.12 104646.77 85219.39 108767.49
Mar 2000       96067.58 88394.02 103741.15 84262.40 107872.77
Apr 2000       81283.94 73589.07  88978.81 69445.98  93121.90
May 2000       72754.64 65037.40  80471.87 60882.27  84627.00
Jun 2000       77627.00 69886.35  85367.64 65718.62  89535.38
Jul 2000       91675.80 83910.71  99440.89 79729.81 103621.80
Aug 2000       83540.16 75749.59  91330.73 71554.97  95525.36
Sep 2000       68577.97 60760.90  76395.04 56552.01  80603.93
Oct 2000       72856.77 65012.19  80701.35 60788.49  84925.05
Nov 2000       83927.02 76053.94  91800.11 71814.89  96039.16
Dec 2000      100091.16 92188.58 107993.74 87933.65 112248.67
accuracy(pronostico)
                       ME     RMSE      MAE        MPE     MAPE      MASE
Training set 5.729326e-13 5050.477 4024.289 -0.3434747 6.726303 0.6524221
                 ACF1
Training set 0.623945
accuracy(pronostico,y.test)
                       ME     RMSE      MAE        MPE     MAPE      MASE
Training set 5.729326e-13 5050.477 4024.289 -0.3434747 6.726303 0.6524221
Test set     4.766925e+03 8704.755 6962.741  4.0412577 7.106520 1.1288072
                  ACF1 Theil's U
Training set 0.6239450        NA
Test set     0.5177485 0.5156388

Extensiones de modelos de regresión

Contenido

  1. Introducción

  2. Modelos de regresión

  3. Regresión espuria

  4. El tiempo como covariables

  5. Pronósticos usando regresión

  6. Extensiones de modelos de regresión

Extensiones de modelos de regresión

  • Modelos no lineales.

  • Modelos lineales generalizados (GLM).

    • Variables dependientes que pertenecen a una familia exponencial: Poisson, Exponencial, etc.
  • Modelos aditivos generalizados (GAM).

    • Efectos fijos y aleatorios.
  • Modelos aditivos generalizados para locación, escala y forma (GAMLSS).

Ejemplo

Temperatura diaria (F) en Cairo de 01-01-1995 hasta 21-05-2005

data(cairo)
plot(cairo$temp~cairo$time,type="l", ylab="temp",xlab="time")

  • La estacionalidad es de 365 o 366 días.

  • Muy costoso crear tantas variables indicadoras.

  • Ajuste de un GAM
ctamm <- gamm(temp~s(day.of.year,bs="cc",k=20)+s(time,bs="cr"),
              data=cairo,correlation=corAR1(form=~1|year))
pred<-fitted(ctamm$gam)
plot(cairo$temp~cairo$time,type="l", ylab="temp",xlab="time")
points(cairo$time,pred,type="l",col=2,lwd=2)
  • Periodograma (análisis espectral)
library(astsa)
x<-cairo$temp
x.spec <-mvspec(x,log="no")

plot(x.spec)
frecuencia<-x.spec$freq[x.spec$spec==max(x.spec$spec)]
abline(v=frecuencia,col=2)

1/frecuencia
[1] 349.0909

La frecuencia dominante \(\omega\) es 0.002864583, es decir, se completa un ciclo aproximadamente cada 349.09 días.

  • Modelo de regresión lineal con covariables senos y cosenos \[y_t=\alpha+\beta_1 cos(2 \pi \omega t)+\beta_2 sen(2 \pi \omega t)+\varepsilon_t\]
cos1<-cos(2*pi*frecuencia*cairo$time)
sin1<-sin(2*pi*frecuencia*cairo$time)
data=data.frame(x=x,time=cairo$time,cos=cos1,sen=sin1)
mod<-lm(x~cos1+sin1,data=data)
plot(cairo$temp~cairo$time,type="l", ylab="temp",xlab="time")
points(cairo$time,fitted(mod),type="l",col=2,lwd=2)

  • Modelo de regresión lineal con covariables senos y cosenos con 5 frecuencias más importantes:

\[y_t=\alpha+ \sum_{k=1}^{5} \left[ \beta_{1k} cos(2 \pi \omega_k t)+\beta_{2k} sen(2 \pi \omega_k t) \right]+\varepsilon_t\]

  • Las 8 frecuencias más importantes:
spectro_orden <- x.spec$details %>% as.data.frame() %>% arrange(desc(spectrum))
spectro_orden[1:8,]
  frequency   period  spectrum
1    0.0029 349.0909 91251.731
2    0.0026 384.0000 54420.640
3    0.0023 426.6667  8986.496
4    0.0031 320.0000  6148.723
5    0.0034 295.3846  2902.172
6    0.0021 480.0000  2409.255
7    0.0018 548.5714  1666.177
8    0.0039 256.0000  1331.423
frecuencia1 <- spectro_orden$frequency[1] 
frecuencia2 <- spectro_orden$frequency[2]
frecuencia3 <- spectro_orden$frequency[3]
frecuencia4 <- spectro_orden$frequency[4]
frecuencia5 <- spectro_orden$frequency[5]
cos1<-cos(2*pi*frecuencia1*cairo$time)
sin1<-sin(2*pi*frecuencia1*cairo$time)
cos2<-cos(2*pi*frecuencia2*cairo$time)
sin2<-sin(2*pi*frecuencia2*cairo$time)
cos3<-cos(2*pi*frecuencia3*cairo$time)
sin3<-sin(2*pi*frecuencia3*cairo$time)
cos4<-cos(2*pi*frecuencia4*cairo$time)
sin4<-sin(2*pi*frecuencia4*cairo$time)
cos5<-cos(2*pi*frecuencia5*cairo$time)
sin5<-sin(2*pi*frecuencia5*cairo$time)

data=data.frame(x=x,
                cos1=cos1,sen1=sin1,
                cos2=cos2,sen2=sin2,
                cos3=cos3,sen3=sin3,
                cos4=cos4,sen4=sin4,
                cos5=cos5,sen5=sin5)

mod2<-lm(x~.,data=data)
plot(cairo$temp~cairo$time,type="l", ylab="temp",xlab="time")
points(cairo$time,fitted(mod2),type="l",col=2,lwd=2)  

Próximo tema

Tema 5: Modelos de series temporales