Actividad ch1

Contexto

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

Las variables son:

• I = Aerolínea
• T = Año
• Q = Millas voladas por los pasajeros
• 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 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)
library(ggplot2)
#install.packages("DataExplorer")
library(DataExplorer)
library(gplots)

Importar la base de datos

df<- read.csv("C:\\Users\\LuisD\\Documents\\OCTAVO SEMESTRE\\Generación de escenarios\\Módulo 1\\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() +
  labs(title =  "Costo por aerolínea (en miles de dólares)", x= "Año", y= "Costo (USD)",
       color= "Aerolínea") + 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= "Índice normalizado",
       color= "Aerolínea") + 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= "Aerolínea") + 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= "Aerolínea") + theme_minimal()

Tema 1. Datos de panel

Heterogeneidad

plotmeans

## function (formula, data = NULL, subset, na.action, bars = TRUE, 
##     p = 0.95, minsd = 0, minbar = NULL, maxbar = NULL, xlab = names(mf)[2], 
##     ylab = names(mf)[1], mean.labels = FALSE, ci.label = FALSE, 
##     n.label = TRUE, text.n.label = "n=", digits = getOption("digits"), 
##     col = "black", barwidth = 1, barcol = "blue", connect = TRUE, 
##     ccol = col, legends = names(means), xaxt, use.t = TRUE, lwd = par("lwd"), 
##     ...) 
## {
##     if (invalid(formula) || (length(formula) != 3)) 
##         stop("formula missing or incorrect")
##     if (invalid(na.action)) 
##         na.action <- options("na.action")
##     mf <- match.call(expand.dots = FALSE)
##     if (is.matrix(eval(mf$data, parent.frame()))) 
##         mf$data <- as.data.frame(data)
##     m <- match(c("formula", "data", "subset", "na.action"), names(mf), 
##         0L)
##     mf <- mf[c(1L, m)]
##     mf[[1L]] <- quote(stats::model.frame)
##     mf <- eval(mf, parent.frame())
##     response <- attr(attr(mf, "terms"), "response")
##     wFact <- which(attr(attr(mf, "terms"), "dataClasses") == 
##         "factor")
##     for (i in wFact) mf[, i] <- factor(mf[, i])
##     means <- sapply(split(mf[[response]], mf[[-response]]), mean, 
##         na.rm = TRUE)
##     ns <- sapply(sapply(split(mf[[response]], mf[[-response]]), 
##         na.omit, simplify = FALSE), length)
##     xlim <- c(0.5, length(means) + 0.5)
##     if (!bars) {
##         plot(means, ..., col = col, xlim = xlim)
##     }
##     else {
##         myvar <- function(x) var(x[!is.na(x)])
##         vars <- sapply(split(mf[[response]], mf[[-response]]), 
##             myvar)
##         vars <- ifelse(vars < (minsd^2), (minsd^2), vars)
##         if (use.t) 
##             ci.width <- qt((1 + p)/2, ns - 1) * sqrt(vars/ns)
##         else ci.width <- qnorm((1 + p)/2) * sqrt(vars/ns)
##         if (length(mean.labels) == 1) {
##             if (mean.labels == TRUE) 
##                 mean.labels <- format(round(means, digits = digits))
##             else if (mean.labels == FALSE) 
##                 mean.labels <- NULL
##         }
##         plotCI(x = 1:length(means), y = means, uiw = ci.width, 
##             xaxt = "n", xlab = xlab, ylab = ylab, labels = mean.labels, 
##             col = col, xlim = xlim, lwd = barwidth, barcol = barcol, 
##             minbar = minbar, maxbar = maxbar, ...)
##         if (invalid(xaxt) || xaxt != "n") 
##             axis(1, at = 1:length(means), labels = legends, ...)
##         if (ci.label) {
##             ci.lower <- means - ci.width
##             ci.upper <- means + ci.width
##             if (!invalid(minbar)) 
##                 ci.lower <- ifelse(ci.lower < minbar, minbar, 
##                   ci.lower)
##             if (!invalid(maxbar)) 
##                 ci.upper <- ifelse(ci.upper > maxbar, maxbar, 
##                   ci.upper)
##             labels.lower <- paste(" \n", format(round(ci.lower, 
##                 digits = digits)), sep = "")
##             labels.upper <- paste(format(round(ci.upper, digits = digits)), 
##                 "\n ", sep = "")
##             text(x = 1:length(means), y = ci.lower, labels = labels.lower, 
##                 col = col)
##             text(x = 1:length(means), y = ci.upper, labels = labels.upper, 
##                 col = col)
##         }
##     }
##     if (n.label) {
##         text(x = 1:length(means), y = par("usr")[3], labels = paste(text.n.label, 
##             ns, "\n", sep = ""))
##     }
##     if (!invalid(connect) & !identical(connect, FALSE)) {
##         if (is.list(connect)) {
##             if (length(ccol) == 1) 
##                 ccol <- rep(ccol, length(connect))
##             for (which in 1:length(connect)) lines(x = connect[[which]], 
##                 y = means[connect[[which]]], col = ccol[which])
##         }
##         else lines(means, ..., lwd = lwd, col = ccol)
##     }
## }
## <bytecode: 0x000001bd62d14520>
## <environment: namespace:gplots>

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

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

Creación de datos panel

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

##              C        Q      PF       LF
## 1-1970 1140640 0.952757  106650 0.534487
## 1-1971 1215690 0.986757  110307 0.532328
## 1-1972 1309570 1.091980  110574 0.547736
## 1-1973 1511530 1.175780  121974 0.540846
## 1-1974 1676730 1.160170  196606 0.591167
## 1-1975 1823740 1.173760  265609 0.575417
## 1-1976 2022890 1.290510  263451 0.594495
## 1-1977 2314760 1.390670  316411 0.597409
## 1-1978 2639160 1.612730  384110 0.638522
## 1-1979 3247620 1.825440  569251 0.676287
## 1-1980 3787750 1.546040  871636 0.605735
## 1-1981 3867750 1.527900  997239 0.614360
## 1-1982 3996020 1.660200  938002 0.633366
## 1-1983 4282880 1.822310  859572 0.650117
## 1-1984 4748320 1.936460  823411 0.625603
## 2-1970  569292 0.520635  103795 0.490851
## 2-1971  640614 0.534627  111477 0.473449
## 2-1972  777655 0.655192  118664 0.503013
## 2-1973  999294 0.791575  114797 0.512501
## 2-1974 1203970 0.842945  215322 0.566782
## 2-1975 1358100 0.852892  281704 0.558133
## 2-1976 1501350 0.922843  304818 0.558799
## 2-1977 1709270 1.000000  348609 0.572070
## 2-1978 2025400 1.198450  374579 0.624763
## 2-1979 2548370 1.340670  544109 0.628706
## 2-1980 3137740 1.326240  853356 0.589150
## 2-1981 3557700 1.248520 1003200 0.532612
## 2-1982 3717740 1.254320  941977 0.526652
## 2-1983 3962370 1.371770  856533 0.540163
## 2-1984 4209390 1.389740  821361 0.528775
## 3-1970  286298 0.262424  118788 0.524334
## 3-1971  309290 0.266433  123798 0.537185
## 3-1972  342056 0.306043  122882 0.582119
## 3-1973  374595 0.325586  131274 0.579489
## 3-1974  450037 0.345706  222037 0.606592
## 3-1975  510412 0.367517  278721 0.607270
## 3-1976  575347 0.409937  306564 0.582425
## 3-1977  669331 0.448023  356073 0.573972
## 3-1978  783799 0.539595  378311 0.654256
## 3-1979  913883 0.539382  555267 0.631055
## 3-1980 1041520 0.467967  850322 0.569240
## 3-1981 1125800 0.450544 1015610 0.589682
## 3-1982 1096070 0.468793  954508 0.587953
## 3-1983 1198930 0.494397  886999 0.565388
## 3-1984 1170470 0.493317  844079 0.577078
## 4-1970  145167 0.086393  114987 0.432066
## 4-1971  170192 0.096740  120501 0.439669
## 4-1972  247506 0.141500  121908 0.488932
## 4-1973  309391 0.169715  127220 0.484181
## 4-1974  354338 0.173805  209405 0.529925
## 4-1975  373941 0.164272  263148 0.532723
## 4-1976  420915 0.170906  316724 0.549067
## 4-1977  474017 0.177840  363598 0.557140
## 4-1978  532590 0.192248  389436 0.611377
## 4-1979  676771 0.242469  547376 0.645319
## 4-1980  880438 0.256505  850418 0.611734
## 4-1981 1052020 0.249657 1011170 0.580884
## 4-1982 1193680 0.273923  951934 0.572047
## 4-1983 1303390 0.371131  881323 0.594570
## 4-1984 1436970 0.421411  831374 0.585525
## 5-1970   91361 0.051028  118222 0.442875
## 5-1971   95428 0.052646  116223 0.462473
## 5-1972   98187 0.056348  115853 0.519118
## 5-1973  115967 0.066953  129372 0.529331
## 5-1974  138382 0.070308  243266 0.557797
## 5-1975  156228 0.073961  277930 0.556181
## 5-1976  183169 0.084946  317273 0.569327
## 5-1977  210212 0.095474  358794 0.583465
## 5-1978  274024 0.119814  397667 0.631818
## 5-1979  356915 0.150046  566672 0.604723
## 5-1980  432344 0.144014  848393 0.587921
## 5-1981  524294 0.169300 1005740 0.616159
## 5-1982  530924 0.172761  958231 0.605868
## 5-1983  581447 0.186670  872924 0.594688
## 5-1984  610257 0.213279  844622 0.635545
## 6-1970   68978 0.037682  117112 0.448539
## 6-1971   74904 0.039784  119420 0.475889
## 6-1972   83829 0.044331  116087 0.500562
## 6-1973   98148 0.050245  122997 0.500344
## 6-1974  118449 0.055046  194309 0.528897
## 6-1975  133161 0.052462  307923 0.495361
## 6-1976  145062 0.056977  323595 0.510342
## 6-1977  170711 0.061490  363081 0.518296
## 6-1978  199775 0.069027  386422 0.546723
## 6-1979  276797 0.092749  564867 0.554276
## 6-1980  381478 0.112640  874818 0.517766
## 6-1981  506969 0.154154 1013170 0.580049
## 6-1982  633388 0.186461  930477 0.556024
## 6-1983  804388 0.246847  851676 0.537791
## 6-1984 1009500 0.304013  819476 0.525775

Modelo 1. Regresión agrupada (Pooled)

# El modelo de regresion agrupada (Pooled) es una técnica de estimacion en datos de panel donde se asume que no hay efectos individuales específicos para cada unidad (Ej. Aerolíneas) ni variaciones en el tiempo, Ignora heterogeneidades
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

# Priebda de Breusch-Pagan (BP): para verificar si el modelo pooled es adecuado.
# Pvalue < 0.05 Avanzamos para usar un modelo de efectos fijos o aleatorio
# Pvalue > 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 Pvalue es > mayor que 0.05 podemos utilizar el modelo pooled.

Modelos 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 4. 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 su R2 ajustadas, el mejor método en el modelo de efectos aleatorios es el 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

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 el Modelo Arima

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_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_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 USD)", xlab= "Año", ylab="USD")

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 análisis factorial confirmatorio

df_escalada <- df
df_escalada$I <- as.numeric(df_escalada$I)
str(df_escalada)

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

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

lavaanPlot(cfa)