Actividad 4

Contexto

La base de datos es de la Universidad de Nueva York y contiene 90 observaciones que incluyen los costos de 6 aerolíneas estdounidenses durante 15 años, de 1970 a 1984.

Las variables son:
* I = Aerolínea
* T = Año
* Q = Millas Voladas por los Pasajeros (Valores Normalizados)
* C = Costo Total en $1,000
* PF = Precio del Combustible
* LF = Factor de Carga (Utilización promedio de la capacidad de la flota)

Fuente: Tabla F7.1

Instalar paquetes y llamar librerías

#install.packages("plm")
library(plm)
#install.packages("tidyverse")
library(tidyverse)
#install.packages("forecast")
library(forecast)
#install.packages("lavaan")
library(lavaan)
#install.packages("lavaanPlot")
library(lavaanPlot)
#install.packages("DataExplorer")
library(DataExplorer)
#install.packages("ggplot2")
library(ggplot2)
#install.packages("gplots")
library(gplots)
library(readxl)

Importar la base de datos

df <- read.csv("C:/Users/ferna/Desktop/Cost Data for U.S. Airlines.csv")

Análisis descriptivo

df$I<- as.factor(df$I)
df$Y<- df$T+1969
summary(df)

##  I            T            C                 Q                 PF         
##  1:15   Min.   : 1   Min.   :  68978   Min.   :0.03768   Min.   : 103795  
##  2:15   1st Qu.: 4   1st Qu.: 292046   1st Qu.:0.14213   1st Qu.: 129848  
##  3:15   Median : 8   Median : 637001   Median :0.30503   Median : 357434  
##  4:15   Mean   : 8   Mean   :1122524   Mean   :0.54499   Mean   : 471683  
##  5:15   3rd Qu.:12   3rd Qu.:1345968   3rd Qu.:0.94528   3rd Qu.: 849840  
##  6:15   Max.   :15   Max.   :4748320   Max.   :1.93646   Max.   :1015610  
##        LF               Y       
##  Min.   :0.4321   Min.   :1970  
##  1st Qu.:0.5288   1st Qu.:1973  
##  Median :0.5661   Median :1977  
##  Mean   :0.5605   Mean   :1977  
##  3rd Qu.:0.5947   3rd Qu.:1981  
##  Max.   :0.6763   Max.   :1984

str(df)

## 'data.frame':    90 obs. of  7 variables:
##  $ I : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ T : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ C : int  1140640 1215690 1309570 1511530 1676730 1823740 2022890 2314760 2639160 3247620 ...
##  $ Q : num  0.953 0.987 1.092 1.176 1.16 ...
##  $ PF: int  106650 110307 110574 121974 196606 265609 263451 316411 384110 569251 ...
##  $ LF: num  0.534 0.532 0.548 0.541 0.591 ...
##  $ Y : num  1970 1971 1972 1973 1974 ...

head(df)

##   I T       C        Q     PF       LF    Y
## 1 1 1 1140640 0.952757 106650 0.534487 1970
## 2 1 2 1215690 0.986757 110307 0.532328 1971
## 3 1 3 1309570 1.091980 110574 0.547736 1972
## 4 1 4 1511530 1.175780 121974 0.540846 1973
## 5 1 5 1676730 1.160170 196606 0.591167 1974
## 6 1 6 1823740 1.173760 265609 0.575417 1975

#create_report(df)
plot_missing(df)

plot_histogram(df)

plot_correlation(df)

ggplot(df,aes(x=Y, y=C, color=I, group=I)) + 
  geom_line() + 
  geom_point()+ 
  labs(title="Costo por Aerolinea (en miles)", x= "Año", y="Costo (USD)", color="Aerolínea") +
  theme_minimal()

ggplot(df,aes(x=Y, y=Q, color=I, group=I)) + 
  geom_line() + 
  geom_point()+ 
  labs(title="Millas voladas por pasajero", x= "Año", y="Indice normalizado", color="Aerolínea") +
  theme_minimal()

ggplot(df,aes(x=Y, y=PF, color=I, group=I)) + 
  geom_line() + 
  geom_point()+ 
  labs(title="Precio del combustible", x= "Año", y="Costo (USD)", color="Aerolínea") +
  theme_minimal()

ggplot(df,aes(x=Y, y=LF, color=I, group=I)) + 
  geom_line() + 
  geom_point()+ 
  labs(title="Factor de carga", x= "Porcentaje", y="Millas", color="Aerolínea") +
  theme_minimal()

plotmeans(C~I,main="Heterogeneidad entre Areolineas", xlab="Aerolínea", ylab="Costo (Miles de USD)", data=df)

## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped

Como el valor promedio (círculo) y el rango intercuartil (lineas azules) varían entre individuos, se observa presencia de heterogeneidad.

Creación de datos de panel

df_panel<-pdata.frame(df, index=c("I","Y"))
df_panel<- df_panel %>% select(-c("I","T","Y"))

Modelo 1. Regresión Agrupada (Pooled)

pooled<-plm(C ~ Q + PF +LF ,data=df_panel,model="pooling")
summary(pooled)

## Pooling Model
## 
## Call:
## plm(formula = C ~ Q + PF + LF, data = df_panel, model = "pooling")
## 
## Balanced Panel: n = 6, T = 15, N = 90
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -520654 -250270   37333  208690  849700 
## 
## Coefficients:
##                Estimate  Std. Error t-value  Pr(>|t|)    
## (Intercept)  1.1586e+06  3.6059e+05  3.2129   0.00185 ** 
## Q            2.0261e+06  6.1807e+04 32.7813 < 2.2e-16 ***
## PF           1.2253e+00  1.0372e-01 11.8138 < 2.2e-16 ***
## LF          -3.0658e+06  6.9633e+05 -4.4027 3.058e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.2647e+14
## Residual Sum of Squares: 6.8177e+12
## R-Squared:      0.94609
## Adj. R-Squared: 0.94421
## F-statistic: 503.118 on 3 and 86 DF, p-value: < 2.22e-16

#Prueba de Breush - Pagan (BP): Para verificar si el modelo es adecuado. 
#p-value <0.05 Avanzamos para usar un Modelo de Efcetos Fijos o Aleatorios

#p-value > 0.05 Podemos usar el modelo Pooled
plmtest(pooled,type="bp")

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  C ~ Q + PF + LF
## chisq = 0.61309, df = 1, p-value = 0.4336
## alternative hypothesis: significant effects

Como el p-value es > 0.05, podemos utilizar el modelo Pooled

Modelo 2. Efectos fijos (within)

within<-plm(C ~ Q + PF +LF ,data=df_panel,model="within")
summary(within)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = C ~ Q + PF + LF, data = df_panel, model = "within")
## 
## Balanced Panel: n = 6, T = 15, N = 90
## 
## Residuals:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -551783 -159259    1796       0  137226  499296 
## 
## Coefficients:
##       Estimate  Std. Error t-value  Pr(>|t|)    
## Q   3.3190e+06  1.7135e+05 19.3694 < 2.2e-16 ***
## PF  7.7307e-01  9.7319e-02  7.9437 9.698e-12 ***
## LF -3.7974e+06  6.1377e+05 -6.1869 2.375e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    5.0776e+13
## Residual Sum of Squares: 3.5865e+12
## R-Squared:      0.92937
## Adj. R-Squared: 0.92239
## F-statistic: 355.254 on 3 and 81 DF, p-value: < 2.22e-16

Modelo 3: Efectos Aleatorios (random) Método Walhus

walhus <- plm(C ~ Q + PF + LF, data=df_panel, model="random", random.method="walhus")
summary(walhus)

## Oneway (individual) effect Random Effect Model 
##    (Wallace-Hussain's transformation)
## 
## Call:
## plm(formula = C ~ Q + PF + LF, data = df_panel, model = "random", 
##     random.method = "walhus")
## 
## Balanced Panel: n = 6, T = 15, N = 90
## 
## Effects:
##                     var   std.dev share
## idiosyncratic 7.339e+10 2.709e+05 0.969
## individual    2.363e+09 4.861e+04 0.031
## theta: 0.1788
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -524180 -243611   39332  199517  824905 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept)  1.1267e+06  3.6994e+05  3.0455  0.002323 ** 
## Q            2.0647e+06  7.1927e+04 28.7051 < 2.2e-16 ***
## PF           1.2075e+00  1.0358e-01 11.6578 < 2.2e-16 ***
## LF          -3.0314e+06  7.1431e+05 -4.2438 2.198e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.0182e+14
## Residual Sum of Squares: 6.5784e+12
## R-Squared:      0.93539
## Adj. R-Squared: 0.93314
## Chisq: 1245.09 on 3 DF, p-value: < 2.22e-16

Modelo 3: Efectos Aleatorios (random) Método Amemiya

amemiya <- plm(C ~ Q + PF + LF, data=df_panel, model="random", random.method="amemiya")
summary(amemiya)

## Oneway (individual) effect Random Effect Model 
##    (Amemiya's transformation)
## 
## Call:
## plm(formula = C ~ Q + PF + LF, data = df_panel, model = "random", 
##     random.method = "amemiya")
## 
## Balanced Panel: n = 6, T = 15, N = 90
## 
## Effects:
##                     var   std.dev share
## idiosyncratic 4.270e+10 2.066e+05 0.084
## individual    4.640e+11 6.812e+05 0.916
## theta: 0.9219
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -603585 -144415   22641  158005  485417 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept)  1.0746e+06  4.2105e+05  2.5522    0.0107 *  
## Q            3.2090e+06  1.6482e+05 19.4695 < 2.2e-16 ***
## PF           8.1014e-01  9.6147e-02  8.4260 < 2.2e-16 ***
## LF          -3.7168e+06  6.1330e+05 -6.0603 1.359e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    5.1238e+13
## Residual Sum of Squares: 3.8227e+12
## R-Squared:      0.92539
## Adj. R-Squared: 0.92279
## Chisq: 1066.71 on 3 DF, p-value: < 2.22e-16

Modelo 3: Efectos Aleatorios (random) Método nerlove

nerlove <- plm(C ~ Q + PF + LF, data=df_panel, model="random", random.method="nerlove")
summary(nerlove)

## Oneway (individual) effect Random Effect Model 
##    (Nerlove's transformation)
## 
## Call:
## plm(formula = C ~ Q + PF + LF, data = df_panel, model = "random", 
##     random.method = "nerlove")
## 
## Balanced Panel: n = 6, T = 15, N = 90
## 
## Effects:
##                     var   std.dev share
## idiosyncratic 3.985e+10 1.996e+05 0.066
## individual    5.602e+11 7.485e+05 0.934
## theta: 0.9313
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -601947 -145039   18713  154903  483623 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept)  1.0752e+06  4.4535e+05  2.4142   0.01577 *  
## Q            3.2323e+06  1.6521e+05 19.5652 < 2.2e-16 ***
## PF           8.0229e-01  9.5804e-02  8.3743 < 2.2e-16 ***
## LF          -3.7338e+06  6.0963e+05 -6.1247 9.084e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    5.1133e+13
## Residual Sum of Squares: 3.7726e+12
## R-Squared:      0.92622
## Adj. R-Squared: 0.92365
## Chisq: 1079.63 on 3 DF, p-value: < 2.22e-16

Comparando sus R2 ajustadas, el mejor método es el modelo de efectos aleatorios de tipo Wallhus

Efectos Fijos vs Efectos Aleatorios

phtest(within,walhus)

## 
##  Hausman Test
## 
## data:  C ~ Q + PF + LF
## chisq = 65.039, df = 3, p-value = 4.919e-14
## alternative hypothesis: one model is inconsistent

Tema 2. Series de Tiempo

Generar serie de tiempo

df_a1<-df[df$I == "1", ]
ts_a1<-ts(df_a1$C, start=1970, frequency = 1)

df_a2<-df[df$I == "2", ]
ts_a2<-ts(df_a2$C, start=1970, frequency = 1)

df_a3<-df[df$I == "3", ]
ts_a3<-ts(df_a3$C, start=1970, frequency = 1)

df_a4<-df[df$I == "4", ]
ts_a4<-ts(df_a4$C, start=1970, frequency = 1)

df_a5<-df[df$I == "5", ]
ts_a5<-ts(df_a5$C, start=1970, frequency = 1)

df_a6<-df[df$I == "6", ]
ts_a6<-ts(df_a6$C, start=1970, frequency = 1)

Generar y comparar Autoarima

arima_a4<- auto.arima(ts_a4)
summary(arima_a4)

## Series: ts_a4 
## ARIMA(0,2,0) 
## 
## sigma^2 = 1.469e+09:  log likelihood = -155.65
## AIC=313.3   AICc=313.66   BIC=313.86
## 
## Training set error measures:
##                    ME    RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 7232.074 35684.5 27472.98 1.761789 5.046326 0.2977402 0.1925091

arima_a1<- auto.arima(ts_a1)
summary(arima_a1)

## Series: ts_a1 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##           drift
##       257691.43
## s.e.   44508.78
## 
## sigma^2 = 2.987e+10:  log likelihood = -188.19
## AIC=380.37   AICc=381.46   BIC=381.65
## 
## Training set error measures:
##                    ME   RMSE      MAE       MPE     MAPE     MASE      ACF1
## Training set 58.86321 160892 129527.1 -1.742419 5.395122 0.502644 0.4084903

arima_a2<- auto.arima(ts_a2)
summary(arima_a2)

## Series: ts_a2 
## ARIMA(0,2,0) 
## 
## sigma^2 = 1.392e+10:  log likelihood = -170.26
## AIC=342.53   AICc=342.89   BIC=343.09
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 11689.89 109830.2 79466.33 1.387268 3.747652 0.3056315 0.3172172

arima_a3<- auto.arima(ts_a3)
summary(arima_a3)

## Series: ts_a3 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##          drift
##       63155.14
## s.e.  13344.11
## 
## sigma^2 = 2.685e+09:  log likelihood = -171.32
## AIC=346.64   AICc=347.74   BIC=347.92
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE     MASE       ACF1
## Training set 14.87618 48235.79 38474.72 -0.9277567 5.324145 0.538349 0.09130379

arima_a5<- auto.arima(ts_a5)
summary(arima_a5)

## Series: ts_a5 
## ARIMA(1,2,0) 
## 
## Coefficients:
##           ar1
##       -0.4543
## s.e.   0.2354
## 
## sigma^2 = 775697764:  log likelihood = -151.09
## AIC=306.18   AICc=307.38   BIC=307.31
## 
## Training set error measures:
##                   ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set 3061.06 24911.01 14171.99 2.393894 4.771228 0.3823654 0.008627682

arima_a6<- auto.arima(ts_a6)
summary(arima_a6)

## Series: ts_a6 
## ARIMA(1,2,0) 
## 
## Coefficients:
##          ar1
##       0.5824
## s.e.  0.2281
## 
## sigma^2 = 386182350:  log likelihood = -146.65
## AIC=297.3   AICc=298.5   BIC=298.43
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set 6829.403 17576.86 10190.16 2.076518 3.550582 0.1516841 -0.2989742

Generar el pronóstico

pronostico_a6<-forecast(arima_a6, level=95,h=5)
pronostico_a6

##      Point Forecast   Lo 95   Hi 95
## 1985        1234478 1195962 1272994
## 1986        1471026 1364365 1577687
## 1987        1714311 1510670 1917953
## 1988        1961521 1635113 2287929
## 1989        2211016 1738872 2683160

plot(pronostico_a6, main="Pronóstico de Costo Total (Miles de doláres) ", xlab="Año", ylab="Dólares")

# Tema 3. Modelo de ecuaciones estructurales

Estructurar el modelo

modelo <- '
          #Regresiones
          C ~ Q + PF + LF + I + Y
          Q ~ PF + I
          PF ~ Y 
          LF ~ I
          # Variables Latentes 
          #Varianzas y Covarianzas 
          C ~~ C 
          Q ~~ Q
          PF ~~ PF 
          LF ~~ LF
          #Intercepto
          '

Generar el análisis Factorial Confirmatorio (CFA)

df_escalada <- df 
df_escalada$I <- as.numeric(df_escalada$I)
df_escalada <- scale(df_escalada)
cfa <- cfa(modelo, df_escalada)
summary(cfa)

## lavaan 0.6-19 ended normally after 3 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        13
## 
##   Number of observations                            90
## 
## Model Test User Model:
##                                                       
##   Test statistic                                63.804
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   C ~                                                 
##     Q                 1.000    0.053   18.826    0.000
##     PF                0.194    0.065    3.000    0.003
##     LF               -0.154    0.025   -6.248    0.000
##     I                 0.105    0.052    1.999    0.046
##     Y                 0.140    0.063    2.211    0.027
##   Q ~                                                 
##     PF                0.239    0.046    5.213    0.000
##     I                -0.871    0.046  -18.985    0.000
##   PF ~                                                
##     Y                 0.931    0.038   24.233    0.000
##   LF ~                                                
##     I                -0.340    0.099   -3.429    0.001
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .C                 0.048    0.007    6.708    0.000
##    .Q                 0.187    0.028    6.708    0.000
##    .PF                0.131    0.020    6.708    0.000
##    .LF                0.875    0.130    6.708    0.000

Generar el modelo de ecuaciones estructurales (SEM)

lavaanPlot(cfa)