Actividad 4. Caso de Negocio 1. Costos en Aerolíneas

#

Contexto

La base de datos es de la universidad de nueva York y contiene 90 observaciones que incluyen los costos de 6 aerolineas estadounidenses durante 15 años, de 1970 a 1984.

Las variables son:

• I = Aerolinea
• T = Año
• Q = Output
• C = Costo total en $1,000
• PF = Precio del Combustible
• LF = Factor de Carga

Instalar paquetes y llamar librerías

library(plm)
library(tidyverse)
library(forecast)
library(lavaan)
library(lavaanPlot)
library(DataExplorer)
library(ggplot2)
library(gplots)

Base de Datos

df <- read.csv("/Users/genarorodriguezalcantara/Desktop/Tec/Generacion de escenarios futuros con analítica (Gpo 101)/PIB/Actividad-4_Caso-de-Negocio-1_Costos-en-Aerolíneas/Cost Data for U.S. Airlines.csv")

Análisis Exploratorio

summary(df)

##        I             T            C                 Q          
##  Min.   :1.0   Min.   : 1   Min.   :  68978   Min.   :0.03768  
##  1st Qu.:2.0   1st Qu.: 4   1st Qu.: 292046   1st Qu.:0.14213  
##  Median :3.5   Median : 8   Median : 637001   Median :0.30503  
##  Mean   :3.5   Mean   : 8   Mean   :1122524   Mean   :0.54499  
##  3rd Qu.:5.0   3rd Qu.:12   3rd Qu.:1345968   3rd Qu.:0.94528  
##  Max.   :6.0   Max.   :15   Max.   :4748320   Max.   :1.93646  
##        PF                LF        
##  Min.   : 103795   Min.   :0.4321  
##  1st Qu.: 129848   1st Qu.:0.5288  
##  Median : 357434   Median :0.5661  
##  Mean   : 471683   Mean   :0.5605  
##  3rd Qu.: 849840   3rd Qu.:0.5947  
##  Max.   :1015610   Max.   :0.6763

str(df)

## 'data.frame':    90 obs. of  6 variables:
##  $ I : int  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 ...

head(df)

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

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

## function (data, output_format = html_document(toc = TRUE, toc_depth = 6, 
##     theme = "yeti"), output_file = "report.html", output_dir = getwd(), 
##     y = NULL, config = configure_report(), report_title = "Data Profiling Report", 
##     ...) 
## {
##     if (!is.data.table(data)) 
##         data <- data.table(data)
##     if (!is.null(y)) {
##         if (!(y %in% names(data))) 
##             stop("`", y, "` not found in data!")
##     }
##     report_dir <- system.file("rmd_template/report.rmd", package = "DataExplorer")
##     suppressWarnings(render(input = report_dir, output_format = output_format, 
##         output_file = output_file, output_dir = output_dir, intermediates_dir = output_dir, 
##         params = list(data = data, report_config = config, response = y, 
##             set_title = report_title), ...))
##     report_path <- path.expand(file.path(output_dir, output_file))
##     browseURL(report_path)
## }
## <bytecode: 0x109977040>
## <environment: namespace:DataExplorer>

Gráficas

#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() +
  labs(title = "Costo por Aerolinea (en miles)", x = "Año", y = "Costo (USD)", color = "Aerolinea") + 
  theme_minimal()

ggplot(df, aes(x = Y, y = Q, color = I, group = I)) + 
  geom_line() +
  labs(title = "Millas Voladas (por pasajero)", x = "Año", y = "Indice Normalizado", color = "Aerolinea") + 
  theme_minimal()

ggplot(df, aes(x = Y, y = PF, color = I, group = I)) + 
  geom_line() +
  labs(title = "Precio del Combustible", x = "Año", y = "Costo (USD)", color = "Aerolinea") + 
  theme_minimal()

ggplot(df, aes(x = Y, y = LF, color = I, group = I)) + 
  geom_line() +
  labs(title = "Factor de Carga", x = "Año", y = "Porcentaje", color = "Aerolinea") + 
  theme_minimal()

# Tema 1. Datos de Panel ## Heterogeneidad

plotmeans(C~I, main = "Heterogeneidad entre Aerolineas", xlab = "Aerolinea", ylab = "Costo (Miles de USD)", data = df)

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

Creacion de Datos de Panel

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

Modelo 1. Regression Agrupada (Pooled)

# El modelo de Regresión Agrupada (pooled) es una técnica de estimación de datos de panel donde se sasume que no hay efectos individuales especificos para cada unidad (Ej. Aerolineas) ni variaciones en el tiempo. Ignora heterogeneidad.
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 Breusch-Pagan (BP): Para verufucar si elmodelo pooled es adecuado.
# p-value < 0.05 Avanzamos para usar un modelo de Efectos 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

Este modelo tiene un p-value mayor a 0.05, lo que nos dice que 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 4. 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 5. 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 en el modelo de efectos aleatorios es el de Walhus.

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

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

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

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

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

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

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

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

## Series: ts_a1 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##           drift
##       257691.43
## s.e.   44509.37
## 
## 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,1,0) with drift 
## 
## Coefficients:
##           drift
##       257691.43
## s.e.   44509.37
## 
## 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_a3 <- auto.arima(ts_a3)
summary(arima_a3)

## Series: ts_a3 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##           drift
##       257691.43
## s.e.   44509.37
## 
## 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_a4 <- auto.arima(ts_a4)
summary(arima_a4)

## Series: ts_a4 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##           drift
##       257691.43
## s.e.   44509.37
## 
## 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_a5 <- auto.arima(ts_a5)
summary(arima_a5)

## Series: ts_a5 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##           drift
##       257691.43
## s.e.   44509.37
## 
## 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_a6 <- auto.arima(ts_a6)
summary(arima_a6)

## Series: ts_a6 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##           drift
##       257691.43
## s.e.   44509.37
## 
## 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

Generar el Pronóstico

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

##      Point Forecast   Lo 95   Hi 95
## 1985        5006011 4667279 5344744
## 1986        5263703 4784663 5742743
## 1987        5521394 4934693 6108096
## 1988        5779086 5101621 6456550
## 1989        6036777 5279349 6794206

plot(pronostico_a6, main="Pronóstico de costo total (en miles)", xlab="Año", ylab="Dólares")

Tema 3. Modelo de Ecuaciones Estructurales

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 Analisis 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 2 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

lavaanPlot(cfa)

Aplicación de Shiny

Link de la Aplicación