Taller 2 Monitorıa: Modelacion de la TendenciaDeterminıstica

1. Problema 1Considere las series en archivo SKYLEPT.DAT y detalladas en la Tabla 1 y Figura 1, y

usando R, para cada una de ellas ajuste un modelo de tendencia conveniente, realice ademasel analisis de residuales.

Tabla 1: Componentes del PIB real, Estados Unidos (millones de $U.S de 1982)Ano AGRICULTURE MANUFACTURING RETAIL SERVICES

1960 68.30 338.70 153.80 190.201961 67.50 339.40 153.20 197.701962 67.10 368.30 163.30 207.701963 67.20 397.40 169.00 217.401964 65.20 425.40 179.40 230.701965 66.70 462.50 190.00 240.401966 62.40 497.90 199.20 253.901967 65.50 496.60 201.50 265.201968 63.60 522.00 211.60 274.701969 65.30 536.70 212.70 287.801970 68.80 506.80 215.60 295.701971 70.60 515.50 224.50 302.401972 70.90 561.20 239.80 320.001973 70.30 621.30 255.60 340.201974 69.70 591.60 245.20 347.501975 73.10 547.50 247.50 352.401976 71.50 600.60 262.80 367.701977 73.30 664.80 270.50 399.601978 73.00 694.70 284.80 421.501979 77.00 712.20 291.30 436.901980 76.40 673.90 281.70 450.901981 87.40 678.60 286.40 463.001982 89.60 634.60 287.50 463.601983 76.70 674.20 307.80 480.401984 84.20 752.40 334.00 509.701985 95.80 779.20 354.40 538.601986 103.60 803.20 377.50 565.801987 105.10 852.20 371.60 592.601988 97.00 917.40 399.20 623.301989 100.50 929.00 412.00 652.30

7080

9010

0

AG

RIC

ULT

UR

E

400

600

800

MA

NU

FAC

TU

RIN

G

150

250

350

RE

TAIL

200

300

400

500

600

1960 1965 1970 1975 1980 1985 1990

SE

RV

ICE

S

Time

datos.stylekp

Figura 1: Series en SKYLEPT.DAT

1

Codigo R 1.1.

#La siguiente linea permite leer los datos guardados en archivo#de texto en disco local. file.choose() abre ventana de windows para#explorar y ubicar el archivo. header=T para indicar que los datos#estan encabezados por nombre de la columna; si no hay encabezado, se toma header=Fdatos.stylekp=read.table(file.choose(),header=T)datos.stylekp=ts(datos.stylekp,frequency=1,start=1960)datos.stylekp

win.graph(width=4.5,height=7,pointsize=8)plot(datos.stylekp)

#Separando las cuatro seriesAGRICULTURE=datos.stylekp[,1]MANUFACTURING=datos.stylekp[,2]RETAIL=datos.stylekp[,3]SERVICES=datos.stylekp[,4]

#definicion de variable ındice de tiempot=1:length(AGRICULTURE);

#O bient=1:nrow(datos.stylekp)

plot(AGRICULTURE)# Ajuste de modelo de tendencia cuadratica para componente AGRICULTURAmodeloA=lm(AGRICULTURE˜t+I(tˆ2))summary(modeloA)

#Creando serie de valores ajustadosyhatA=ts(fitted(modeloA),frequency=1,start=1960)

#Graficando la serie y su ajusteplot(AGRICULTURE)lines(yhatA,col=2,lty=2)

# Grafico de residuales comunes del componente de Agricultura versus variable respuestaplot(fitted(modeloA), residuals(modeloA), type="p", main="Ajustados vs.Residules comunes")abline(h=0)

# Grafico de residuales comunes del componente de Agricultura versus tiempoplot(t, residuals(modeloA), type="l", main="Tiempo vs. Residuales comunes")abline(h=0)

# Prueba de NormalidadtestA=shapiro.test(residuals(modeloA))testA

qqnorm(residuals(modeloA))qqline(residuals(modeloA), main="Prueba de normalidad de residuales de modelo componente Agricultura")legend("topleft", legend=rbind(names(testA),testA), cex=0.8)

plot(MANUFACTURING)# Ajuste de modelo de tendencia cubica para componente MANUFACTURAmodeloM=lm(MANUFACTURING˜t+I(tˆ2)+I(tˆ3))summary(modeloM)

#Creando serie de valores ajustadosyhatM=ts(fitted(modeloM),frequency=1,start=1960)

#Graficando la serie y su ajusteplot(MANUFACTURING)lines(yhatM,col=2,lty=2)

# Grafico de residuales comunes del componente de Manufactura versus variable respuestaplot(fitted(modeloM), residuals(modeloM), type="p", main="Ajustados vs.Residules comunes")abline(h=0)

# Grafico de residuales comunes del componente de Manufactura versus tiempoplot(t, residuals(modeloM), type="l", main="Tiempo vs. Residuales comunes")abline(h=0)

# Prueba de NormalidadtestM=shapiro.test(residuals(modeloM))testM

qqnorm(residuals(modeloM))

2

qqline(residuals(modeloM), main="Prueba de normalidad de residuales de modelo componente Manufactura")legend("topleft", legend=rbind(names(testM),testM), cex=0.8)

plot(RETAIL)# Ajuste de modelo de tendencia cubica para componente RETAILmodeloR=lm(RETAIL˜t+I(tˆ2)+I(tˆ3))summary(modeloR)

#Creando serie de valores ajustadosyhatR=ts(fitted(modeloR),frequency=1,start=1960)

#Graficando la serie y su ajusteplot(RETAIL)lines(yhatR,col=2,lty=2)

# Grafico de residuales comunes del componente de Retail versus variable respuestaplot(fitted(modeloR), residuals(modeloR), type="p", main="Ajustados vs.Residules comunes")abline(h=0)

# Grafico de residuales comunes del componente de Retail versus tiempoplot(t,residuals(modeloR), type="l", main="Tiempo vs. Residuales comunes")abline(h=0)

# Prueba de NormalidadtestR=shapiro.test(rstudent(modeloR))testR

qqnorm(residuals(modeloR))qqline(residuals(modeloR), main="Prueba de normalidad de residuales de modelo componente Retail")legend("topleft", legend=rbind(names(testR),testR), cex=0.8)

plot(SERVICES)# Ajuste de modelo de tendencia cubica para componente SERVICIOSmodeloS=lm(SERVICES˜t+I(tˆ2)+I(tˆ3))summary(modeloS)

#Creando serie de valores ajustadosyhatS=ts(fitted(modeloS),frequency=1,start=1960)

#Graficando la serie y su ajusteplot(SERVICES)lines(yhatS,col=2,lty=2)

# Grafico de residuales comunes del componente de Servicios versus variable respuestaplot(fitted(modeloS), residuals(modeloS), type="p", main="Ajustados vs.Residules comunes")abline(h=0)

# Grafico de residuales comunes del componente de Servicios versus tiempoplot(t, residuals(modeloS), type="l", main="Tiempo vs. Residuales comunes")abline(h=0)

# Prueba de NormalidadtestS=shapiro.test(residuals(modeloS))testS

qqnorm(residuals(modeloS))qqline(residuals(modeloS), main="Prueba de normalidad de residuales de modelo componente Servicios")legend("topleft", legend=rbind(names(testS),testS), cex=0.8)

#Graficando densidad de los residualesplot(density(residuals(modeloS)))

2. Problema 2Los datos que se muestran en la Tabla 2 y Figura 2 pertenecen al ındice de productividad

mensual en Canada, desde Enero de 1950 hasta Diciembre de 1973, y estan guardados enarchivo IndProductividad.txt.

3

Tabla 2: Indice de productividad mensual en CanadaAno ene feb mar abr may jun jul ago sep oct nov dic

1950 77.5 77.6 78.1 78.3 78.3 78.9 79.5 80.0 80.7 82.0 82.4 82.51951 83.4 84.4 85.8 86.5 86.8 88.0 88.7 89.4 90.2 90.6 91.3 91.41952 91.5 91.0 90.5 90.4 89.7 89.8 89.9 89.8 89.9 89.8 89.9 89.61953 89.6 89.4 88.9 88.7 88.5 88.9 89.3 89.6 89.9 90.3 89.9 89.61954 89.6 89.6 89.4 89.5 89.4 89.9 89.9 90.6 90.4 90.4 90.4 90.21955 90.1 90.0 89.8 89.9 90.1 89.7 89.8 90.1 90.4 90.5 90.5 90.51956 90.4 90.1 90.1 90.2 90.2 91.2 91.7 92.2 92.1 92.7 93.1 93.21957 93.1 93.3 93.3 93.6 93.7 94.1 94.3 94.9 95.4 95.5 95.4 95.31958 95.5 95.7 96.2 96.9 96.8 96.8 96.5 96.9 97.2 97.5 97.8 97.71959 97.6 97.3 97.1 97.1 97.2 96.7 97.4 97.8 98.4 99.1 99.3 99.01960 98.7 98.5 98.2 98.7 98.6 98.8 98.7 99.0 99.4 100.2 100.3 100.31961 100.0 99.8 99.9 99.9 99.8 99.8 99.8 99.9 99.9 100.0 100.4 100.51962 100.4 100.5 100.4 100.9 100.7 101.0 101.4 101.7 101.4 101.8 102.1 102.11963 102.2 102.2 102.2 102.4 102.4 102.8 103.3 103.6 103.3 103.4 103.7 103.91964 103.9 104.1 104.2 104.5 104.5 104.7 105.4 105.3 105.0 105.0 105.2 105.91965 106.0 106.2 106.3 106.6 106.8 107.6 108.0 107.9 107.7 107.8 108.5 109.01966 109.3 110.0 110.2 110.8 111.0 111.3 111.7 112.2 112.3 112.5 112.6 112.91967 113.0 113.1 113.4 114.4 114.6 115.2 116.3 116.8 116.6 116.5 116.9 117.51968 118.1 118.2 118.6 119.3 119.3 119.7 120.4 120.7 121.1 121.4 121.9 122.31969 122.6 122.6 123.2 124.6 124.9 125.9 126.4 126.9 126.6 126.8 127.4 127.91970 128.2 128.7 128.9 129.7 129.6 129.9 130.5 130.5 130.2 130.3 130.3 129.81971 130.3 130.9 131.3 132.2 132.7 133.0 134.1 135.0 134.7 134.9 135.4 136.31972 136.7 137.3 137.4 138.2 138.3 138.5 140.2 141.3 141.8 142.0 142.3 143.31973 144.5 145.3 145.7 147.3 148.4 149.7 151.0 153.0 153.9 154.3 155.5 156.4

Time

dato

s.in

dpro

1950 1955 1960 1965 1970

8010

012

014

0

Figura 2: Serie en IndProductividad.txt

Del analisis grafico de la serie, responda

1. ¿Que tipo de tendencia tiene esta serie?

2. Proponga modelos para la serie completa y ajustelos usando R. Concluya si estos mo-delos son apropiados a partir de los graficos de residuales, y sus medidas de ajuste.En cada caso, escriba la ecuacion y el modelo ajustado. Evalue tambien los graficos deresiduales en terminos de la validez del supuesto de independencia.

3. Considere ahora la serie sin sus primeros cinco anos, es decir, desde enero de 1955, yrepita 2.

4. ¿Que se concluye de 2) y 3) en terminos del ajuste de los modelos y comportamientode la serie? ¿las observaciones de los primeros cinco anos afectan significativamente lamodelacion de la serie? ¿por que?

4

Codigo R 2.1.

#La siguiente linea permite leer los datos guardados en archivo#de texto en disco local. file.choose() abre ventana de windows para#explorar y ubicar el archivo. header=T para indicar que los datos#estan encabezados por nombre de la columna; si no hay encabezado, se toma header=F

datos.indpro=scan(file.choose(),dec = ".")datos.indpro=ts(datos.indpro,frequency=12,start=c(1950,1))datos.indpro

#Analisis descriptivowin.graph()plot(datos.indpro)

plot(decompose(datos.indpro,type="additive"))

win.graph(width=5,height=4,pointsize=8)boxplot(datos.indpro˜cycle(datos.indpro),names=month.abb)

#Creando ındice t de tiempot=1:length(datos.indpro)

#Ajuste modelo de tendencia cuadraticamodelo1=lm(datos.indpro˜t+I(tˆ2))summary(modelo1)ajuste=ts(fitted(modelo1),frequency=12,start=c(1950,1))

#Graficando la serie real y la ajustada en mismo planoplot(datos.indpro)lines(ajuste,col=2)legend("topleft",legend=c("Serie real","Serie Ajustada"),col=c(1,2),lty=1)

# Grafico de residuales comunes versus variable respuestaplot(fitted(modelo1), residuals(modelo1), type="p", main="Ajustados vs.Residules comunes")abline(h=0,col=2)

# Grafico de residuales comunes versus tiempoplot(t, residuals(modelo1), type="l", main="Tiempo vs. Residuales comunes")abline(h=0,col=2)

# Prueba de Normalidadtest1=shapiro.test(residuals(modelo1))test1

qqnorm(residuals(modelo1))qqline(residuals(modelo1), main="Prueba de normalidad de residuales de modelo tendencia cuadratica")legend("topleft", legend=rbind(names(test1),test1), cex=0.8)

#Ajuste modelo de tendencia cubicamodelo2=lm(datos.indpro˜t+I(tˆ2)+I(tˆ3))summary(modelo2)ajuste2=ts(fitted(modelo2),frequency=12,start=c(1950,1))

#Graficando la serie real y la ajustada en mismo planoplot(datos.indpro)lines(ajuste2,col=2)legend("topleft",legend=c("Serie real","Serie Ajustada"),col=c(1,2),lty=1)

# Grafico de residuales comunes versus variable respuestaplot(fitted(modelo2), residuals(modelo2), type="p", main="Ajustados vs.Residules comunes")abline(h=0,col=2)

# Grafico de residuales comunes versus tiempoplot(t, residuals(modelo2), type="l", main="Tiempo vs. Residuales comunes")abline(h=0,col=2)

# Prueba de Normalidadtest2=shapiro.test(residuals(modelo2))test2

qqnorm(residuals(modelo2))qqline(residuals(modelo2), main="Prueba de normalidad de residuales de modelo tendencia exponencial")legend("topleft", legend=rbind(names(test2),test2), cex=0.8)

#Ajuste modelo de tendencia polinomio p=4modelo3=lm(datos.indpro˜t+I(tˆ2)+I(tˆ3)+I(tˆ4))summary(modelo3)

5

ajuste3=ts(fitted(modelo3),frequency=12,start=c(1950,1))plot(datos.indpro)lines(ajuste3,col=2)legend("topleft",legend=c("Serie real","Serie Ajustada"),col=c(1,2),lty=1)

AIC(modelo1)AIC(modelo2)AIC(modelo3)BIC(modelo1)BIC(modelo2)BIC(modelo3)

#Parte 2datos.indpro2=ts(datos.indpro[61:length(datos.indpro)],freq=12,start=c(1955,1))

plot(datos.indpro2)

#Creando ındice t de tiempot=1:length(datos.indpro2)

#Ajuste modelo de tendencia cuadraticamodelo1b=lm(datos.indpro2˜t+I(tˆ2))summary(modelo1b)ajusteb=ts(fitted(modelo1b),frequency=12,start=c(1955,1))

plot(datos.indpro2)lines(ajusteb,col=2)legend("topleft",legend=c("Serie real","Serie Ajustada"),col=c(1,2),lty=1)

# Grafico de residuales comunes versus variable respuestaplot(fitted(modelo1b), residuals(modelo1b), type="p", main="Ajustados vs.Residules comunes")abline(h=0,col=2)

# Grafico de residuales comunes versus tiempoplot(t, residuals(modelo1b), type="l", main="Tiempo vs. Residuales comunes")abline(h=0,col=2)

# Prueba de Normalidadtest1b=shapiro.test(residuals(modelo1b))test1b

qqnorm(residuals(modelo1b))qqline(residuals(modelo1b), main="Prueba de normalidad de residuales de modelo tendencia cuadratica")legend("topleft", legend=rbind(names(test1b),test1b), cex=0.8)

#Correr modelo cubico para estos datos y comparar con el modelo cuadratico con AIC, BIC y comportamiento de residuosAIC(modelo1b)BIC(modelo1b)

6