rm(list=ls())
#setwd("C:/DAVE2/CIDE/Doc/Econometria II/Lab_8B")
options(scipen=999) 

library(AER)  #Econometría Aplicada con R

## Loading required package: car

## Loading required package: carData

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: survival

library(plm)  #Ajuste de modelos con datos panel
library(stargazer) # Para salidas de las regresiones

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

library(tseries) # for `adf.test()`

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(dynlm) #for function `dynlm()`
library(vars) # for function `VAR()`

## Loading required package: MASS

## Loading required package: strucchange

## Loading required package: urca

library(nlWaldTest) # for the `nlWaldtest()` function
library(lmtest) #for `coeftest()` and `bptest()`.
library(broom) #for `glance(`) and `tidy()`
library(car) #for `hccm()` robust standard errors
library(sandwich)
library(knitr) #for `kable()`
library(forecast) 
library(systemfit)

## Loading required package: Matrix

## 
## Please cite the 'systemfit' package as:
## Arne Henningsen and Jeff D. Hamann (2007). systemfit: A Package for Estimating Systems of Simultaneous Equations in R. Journal of Statistical Software 23(4), 1-40. http://www.jstatsoft.org/v23/i04/.
## 
## If you have questions, suggestions, or comments regarding the 'systemfit' package, please use a forum or 'tracker' at systemfit's R-Forge site:
## https://r-forge.r-project.org/projects/systemfit/

library(AER)
library(xtable)
library(readxl)
library(ggplot2)

Manejo de datos de Panel

Problema

Queremos convertir nuestros datos de wide a long

Muchas funciones en R pueden hacer eso. Lo necesitamos a veces para modelos de panel. series de tiempo y modelos multinomiales.

Cuidado con el formato necesario para cada paquete.

Podemos utilizar dos librerias distintas:

gather() y spread() de tydr
melt() y dcast() del paquete reshape

Creaciòn de datos

Estos datasets contienen los mismos datos, pero en formatos amplios y largos. Cada uno de ellos se convertirá al otro formato a continuación.

olddata_wide <- read.table(header=TRUE, text='
 subject sex control cond1 cond2
       1   M     7.9  12.3  10.7
       2   F     6.3  10.6  11.1
       3   F     9.5  13.1  13.8
       4   M    11.5  13.4  12.9
')
# Make sure the subject column is a factor
olddata_wide$subject <- factor(olddata_wide$subject)

Y ahora la tabla long:

olddata_long <- read.table(header=TRUE, text='
 subject sex condition measurement
       1   M   control         7.9
       1   M     cond1        12.3
       1   M     cond2        10.7
       2   F   control         6.3
       2   F     cond1        10.6
       2   F     cond2        11.1
       3   F   control         9.5
       3   F     cond1        13.1
       3   F     cond2        13.8
       4   M   control        11.5
       4   M     cond1        13.4
       4   M     cond2        12.9
')
# Make sure the subject column is a factor
olddata_long$subject <- factor(olddata_long$subject)

Empezamos con tydr y con la transformacion: wide -> long

Utilizando gather:

library(tidyr)

## 
## Attaching package: 'tidyr'

## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack

Ploteamos:

olddata_wide

##   subject sex control cond1 cond2
## 1       1   M     7.9  12.3  10.7
## 2       2   F     6.3  10.6  11.1
## 3       3   F     9.5  13.1  13.8
## 4       4   M    11.5  13.4  12.9

Transformamos:

# The arguments to gather():
# - data: Data object
# - key: Name of new key column (made from names of data columns)
# - value: Name of new value column
# - ...: Names of source columns that contain values
# - factor_key: Treat the new key column as a factor (instead of character vector)
data_long <- gather(olddata_wide, condition, measurement, control:cond2, factor_key=TRUE)
data_long

##    subject sex condition measurement
## 1        1   M   control         7.9
## 2        2   F   control         6.3
## 3        3   F   control         9.5
## 4        4   M   control        11.5
## 5        1   M     cond1        12.3
## 6        2   F     cond1        10.6
## 7        3   F     cond1        13.1
## 8        4   M     cond1        13.4
## 9        1   M     cond2        10.7
## 10       2   F     cond2        11.1
## 11       3   F     cond2        13.8
## 12       4   M     cond2        12.9

En este ejemplo, las columnas de origen que se recopilan se especifican con control: cond2. Esto significa utilizar todas las columnas, posicionalmente, entre control y cond2. Otra forma de hacerlo es nombrar las columnas individualmente, como en:

gather(olddata_wide, condition, measurement, control, cond1, cond2)

Si necesita utilizar gather() mediante programación, es posible que deba utilizar variables que contengan nombres de columna. Para hacer esto, debe usar la función gather_() en su lugar, que toma cadenas en lugar de nombres de columnas (sin comillas).

keycol <- "condition"
valuecol <- "measurement"
gathercols <- c("control", "cond1", "cond2")

gather_(olddata_wide, keycol, valuecol, gathercols)

##    subject sex condition measurement
## 1        1   M   control         7.9
## 2        2   F   control         6.3
## 3        3   F   control         9.5
## 4        4   M   control        11.5
## 5        1   M     cond1        12.3
## 6        2   F     cond1        10.6
## 7        3   F     cond1        13.1
## 8        4   M     cond1        13.4
## 9        1   M     cond2        10.7
## 10       2   F     cond2        11.1
## 11       3   F     cond2        13.8
## 12       4   M     cond2        12.9

Utilizamos ahroa spread para pasar de long -> wide

olddata_long

##    subject sex condition measurement
## 1        1   M   control         7.9
## 2        1   M     cond1        12.3
## 3        1   M     cond2        10.7
## 4        2   F   control         6.3
## 5        2   F     cond1        10.6
## 6        2   F     cond2        11.1
## 7        3   F   control         9.5
## 8        3   F     cond1        13.1
## 9        3   F     cond2        13.8
## 10       4   M   control        11.5
## 11       4   M     cond1        13.4
## 12       4   M     cond2        12.9

Ahora transformamos:

# The arguments to spread():
# - data: Data object
# - key: Name of column containing the new column names
# - value: Name of column containing values
data_wide <- spread(olddata_long, condition, measurement)
data_wide

##   subject sex cond1 cond2 control
## 1       1   M  12.3  10.7     7.9
## 2       2   F  10.6  11.1     6.3
## 3       3   F  13.1  13.8     9.5
## 4       4   M  13.4  12.9    11.5

Para obtener un dataset màsa bonito podemos:

# Rename cond1 to first, and cond2 to second
names(data_wide)[names(data_wide)=="cond1"] <- "first"
names(data_wide)[names(data_wide)=="cond2"] <- "second"

# Reorder the columns
data_wide <- data_wide[, c(1,2,5,3,4)]
data_wide

##   subject sex control first second
## 1       1   M     7.9  12.3   10.7
## 2       2   F     6.3  10.6   11.1
## 3       3   F     9.5  13.1   13.8
## 4       4   M    11.5  13.4   12.9

#>   subject sex control first second
#> 1       1   M     7.9  12.3   10.7
#> 2       2   F     6.3  10.6   11.1
#> 3       3   F     9.5  13.1   13.8
#> 4       4   M    11.5  13.4   12.9

Ahora pdoemos utilizar tambien:

df %>% pivot_wider(names_from = key, values_from = value)

Con el paquete reshape2

Cuando queremos pasar de wide -> long

olddata_wide

##   subject sex control cond1 cond2
## 1       1   M     7.9  12.3  10.7
## 2       2   F     6.3  10.6  11.1
## 3       3   F     9.5  13.1  13.8
## 4       4   M    11.5  13.4  12.9

Con nuestra nueva librearia:

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

# Specify id.vars: the variables to keep but not split apart on
melt(olddata_wide, id.vars=c("subject", "sex"))

##    subject sex variable value
## 1        1   M  control   7.9
## 2        2   F  control   6.3
## 3        3   F  control   9.5
## 4        4   M  control  11.5
## 5        1   M    cond1  12.3
## 6        2   F    cond1  10.6
## 7        3   F    cond1  13.1
## 8        4   M    cond1  13.4
## 9        1   M    cond2  10.7
## 10       2   F    cond2  11.1
## 11       3   F    cond2  13.8
## 12       4   M    cond2  12.9

Hay una opcion de melt que permite trabajar uin poco mas facilmente con los outputs:

data_long <- melt(olddata_wide,
        # ID variables - all the variables to keep but not split apart on
    id.vars=c("subject", "sex"),
        # The source columns
    measure.vars=c("control", "cond1", "cond2" ),
        # Name of the destination column that will identify the original
        # column that the measurement came from
    variable.name="condition",
    value.name="measurement"
)
data_long

##    subject sex condition measurement
## 1        1   M   control         7.9
## 2        2   F   control         6.3
## 3        3   F   control         9.5
## 4        4   M   control        11.5
## 5        1   M     cond1        12.3
## 6        2   F     cond1        10.6
## 7        3   F     cond1        13.1
## 8        4   M     cond1        13.4
## 9        1   M     cond2        10.7
## 10       2   F     cond2        11.1
## 11       3   F     cond2        13.8
## 12       4   M     cond2        12.9

Si omitimos measure.vars, melt utilizará automáticamente todas las otras variables como id.vars. Lo contrario es cierto si omite id.vars.

Si no especifica variable.name, nombrará esa columna como “variable”, y si omite value.name, nombrará esa columna como “measurement”.

Opcional: podemos cambiar el nombre de los niveles de factor de la columna de variable.

# Rename factor names from "cond1" and "cond2" to "first" and "second"
levels(data_long$condition)[levels(data_long$condition)=="cond1"] <- "first"
levels(data_long$condition)[levels(data_long$condition)=="cond2"] <- "second"

# Sort by subject first, then by condition
data_long <- data_long[ order(data_long$subject, data_long$condition), ]
data_long

##    subject sex condition measurement
## 1        1   M   control         7.9
## 5        1   M     first        12.3
## 9        1   M    second        10.7
## 2        2   F   control         6.3
## 6        2   F     first        10.6
## 10       2   F    second        11.1
## 3        3   F   control         9.5
## 7        3   F     first        13.1
## 11       3   F    second        13.8
## 4        4   M   control        11.5
## 8        4   M     first        13.4
## 12       4   M    second        12.9

Ahora de long -> wide utilizando dcast()

olddata_long

##    subject sex condition measurement
## 1        1   M   control         7.9
## 2        1   M     cond1        12.3
## 3        1   M     cond2        10.7
## 4        2   F   control         6.3
## 5        2   F     cond1        10.6
## 6        2   F     cond2        11.1
## 7        3   F   control         9.5
## 8        3   F     cond1        13.1
## 9        3   F     cond2        13.8
## 10       4   M   control        11.5
## 11       4   M     cond1        13.4
## 12       4   M     cond2        12.9

Transformamos:

# From the source:
# "subject" and "sex" are columns we want to keep the same
# "condition" is the column that contains the names of the new column to put things in
# "measurement" holds the measurements
library(reshape2)

data_wide <- dcast(olddata_long, subject + sex ~ condition, value.var="measurement")
data_wide

##   subject sex cond1 cond2 control
## 1       1   M  12.3  10.7     7.9
## 2       2   F  10.6  11.1     6.3
## 3       3   F  13.1  13.8     9.5
## 4       4   M  13.4  12.9    11.5

Para que los datos se vean mejor:

# Rename cond1 to first, and cond2 to second
names(data_wide)[names(data_wide)=="cond1"] <- "first"
names(data_wide)[names(data_wide)=="cond2"] <- "second"

# Reorder the columns
data_wide <- data_wide[, c(1,2,5,3,4)]
data_wide

##   subject sex control first second
## 1       1   M     7.9  12.3   10.7
## 2       2   F     6.3  10.6   11.1
## 3       3   F     9.5  13.1   13.8
## 4       4   M    11.5  13.4   12.9

Ejemplo 1

 install.packages("devtools")  # if not already installed
 library(devtools)
 install_git("https://github.com/ccolonescu/PoEdata")

Los datos de panel recopilan información sobre varios individuos (unidades transversales) durante varios períodos.
El panel está equilibrado si se observan todas las unidades en todos los períodos; si faltan algunas unidades en algunos períodos, el panel está desequilibrado.
La ecuación 1 da la forma de un modelo de datos de panel agrupado, donde el subíndice i = 1, …, N denota un individuo (unidad de sección transversal), y t = 1, …, T denota el período de tiempo, o longitudinal unidad. El número total de observaciones en el panel es N × T.

$$ y_{it}={1}+{2}x_{2it}+…+{K}x{Kit}+e_{it}…(1)

$$ ### Los datos

Un panel ancho tiene la dimensión de la sección transversal (N) mucho mayor que la dimensión longitudinal (T); cuando ocurre lo contrario, tenemos un panel largo. Normalmente, se observan las mismas unidades en todos los períodos; cuando este no es el caso y cada período muestrea principalmente otras unidades, el resultado no es un panel de datos adecuado, sino un modelo de secciones transversales agrupadas.

Utilizaremos el paquete plm que nos da tambien la opcion de organizar los datos con forma de pane.

Los datos pueden ser organizado de forma long o wide:

La forma larga tiene una columna para cada variable y una fila para cada período individual; la forma ancha tiene una columna para cada período variable y una fila para cada individuo. La mayoría de los métodos de datos de panel requieren el formato largo, pero muchas fuentes de datos proporcionan una tabla de formato amplio para cada variable.

La siguiente secuencia de código crea una estructura de panel para el conjunto de datos nls_panel utilizando la función pdata.frame del paquete plm y muestra una pequeña parte de este conjunto de datos.

Cargamos la libreria:

library("PoEdata")

Cargamos los datos:

data("nls_panel", package="PoEdata")
nlspd <- pdata.frame(nls_panel, index=c("id", "year"))
smpl <- nlspd[nlspd$id %in% c(1,2),c(1:6, 14:15)]
tbl <- xtable(smpl) 
kable(tbl, digits=4, align="c",
      caption="A data sample")

A data sample
	id	year	lwage	hours	age	educ	union	exper
1-82	1	82	1.8083	38	30	12	1	7.6667
1-83	1	83	1.8634	38	31	12	1	8.5833
1-85	1	85	1.7894	38	33	12	1	10.1795
1-87	1	87	1.8465	40	35	12	1	12.1795
1-88	1	88	1.8564	40	37	12	1	13.6218
2-82	2	82	1.2809	48	36	17	0	7.5769
2-83	2	83	1.5159	43	37	17	0	8.3846
2-85	2	85	1.9302	35	39	17	0	10.3846
2-87	2	87	1.9190	42	41	17	1	12.0385
2-88	2	88	2.2010	42	43	17	1	13.2115

la funciòn pdim() extrae las dimension del panel:

pdim(nlspd)

## Balanced Panel: n = 716, T = 5, N = 3580

El modelo pooled

Un modelo agrupado tiene la especificación de la Ecuación 1, que no permite diferencias de intercepto o pendiente entre individuos. Dicho modelo se puede estimar en R utilizando la combinación de especificaciones en la función plm (), como ilustra la siguiente secuencia de código.

wage.pooled <- plm(lwage~educ+exper+I(exper^2)+
  tenure+I(tenure^2)+black+south+union, 
  model="pooling", data=nlspd)
kable(tidy(wage.pooled), digits=3, 
           caption="Pooled model")

Pooled model
term	estimate	std.error	statistic	p.value
(Intercept)	0.477	0.056	8.487	0.000
educ	0.071	0.003	26.567	0.000
exper	0.056	0.009	6.470	0.000
I(exper^2)	-0.001	0.000	-3.176	0.002
tenure	0.015	0.004	3.394	0.001
I(tenure^2)	0.000	0.000	-1.886	0.059
black	-0.117	0.016	-7.426	0.000
south	-0.106	0.014	-7.465	0.000
union	0.132	0.015	8.839	0.000

La función plm () acepta los siguientes argumentos principales, donde los parámetros mostrados como vectores c (…), como efecto y modelo, solo pueden tomar un valor a la vez de la lista proporcionada.

plm(formula, data, subset, na.action, effect = c("individual", "time", "twoways"), model = c("within", "random", "ht", "between", "pooling", "fd"),...)

Si queremos errores robustos y clusterizados:

tbl <- tidy(coeftest(wage.pooled, vcov=vcovHC(wage.pooled,
                    type="HC0",cluster="group")))
kable(tbl, digits=5, caption=
"Pooled 'wage' model with cluster robust standard errors")

Pooled ‘wage’ model with cluster robust standard errors
term	estimate	std.error	statistic	p.value
(Intercept)	0.47660	0.08441	5.64629	0.00000
educ	0.07145	0.00549	13.01550	0.00000
exper	0.05569	0.01129	4.93242	0.00000
I(exper^2)	-0.00115	0.00049	-2.33440	0.01963
tenure	0.01496	0.00711	2.10401	0.03545
I(tenure^2)	-0.00049	0.00041	-1.18697	0.23532
black	-0.11671	0.02808	-4.15602	0.00003
south	-0.10600	0.02701	-3.92422	0.00009
union	0.13224	0.02703	4.89327	0.00000

Efectos fijos

El modelo de efectos fijos tiene en cuenta las diferencias individuales, traducidas en diferentes intersecciones de la línea de regresión para diferentes individuos. En este caso, el modelo asigna el subíndice i al término constante $\beta_1$ como se muestra en la Ecuación 2; los términos constantes calculados de esta manera se denominan efectos fijos.

\[y_{it}=\beta_{1i}+\beta_{2i}x_{2it}+\beta_{3i}x_{3it}+e_{it}...(2)\] Las variables que cambian poco o nada a lo largo del tiempo, como algunas características individuales, no deben incluirse en un modelo de efectos fijos porque producen colinealidad con los efectos fijos.

nls10 <- pdata.frame(nls_panel[nls_panel$id %in% 1:10,])
## series nev_mar, not_smsa, south are constants and have been removed
wage.fixed <- lm(lwage~exper+I(exper^2)+
                  tenure+I(tenure^2)+union+factor(id)-1,
                  data=nls10)
kable(tidy(wage.fixed), digits=3, 
      caption="Fixed effects in a subsample")

Fixed effects in a subsample
term	estimate	std.error	statistic	p.value
exper	0.238	0.188	1.268	0.213
I(exper^2)	-0.008	0.008	-1.036	0.307
tenure	-0.012	0.034	-0.362	0.720
I(tenure^2)	0.002	0.003	0.854	0.399
union	0.114	0.151	0.753	0.457
factor(id)1	0.152	1.097	0.139	0.891
factor(id)2	0.187	1.071	0.174	0.863
factor(id)3	-0.063	1.351	-0.047	0.963
factor(id)4	0.186	1.343	0.138	0.891
factor(id)5	0.939	1.098	0.855	0.398
factor(id)6	0.794	1.112	0.715	0.480
factor(id)7	0.581	1.236	0.470	0.641
factor(id)8	0.538	1.097	0.490	0.627
factor(id)9	0.418	1.084	0.386	0.702
factor(id)10	0.615	1.090	0.564	0.577

La Tabla precedente muestra los resultados de una regresión OLS en una submuestra de los primeros 10 individuos en el conjunto de datos nls_panel.

La novedad es usar la variable de factor id. La función factor () genera variables ficticias para todas las categorías de la variable, tomando la primera categoría como referencia. Para incluir la referencia en la salida, es necesario excluir la constante del modelo de regresión incluyendo el término -1 en la fórmula de regresión. Cuando no se excluye la constante, los coeficientes de las variables ficticias representan, como es habitual, la diferencia entre la categoría respectiva y la de referencia.

Sin embargo, para estimar un efecto fijo en R no es necesario crear las variables ficticias, sino usar la opción model = “within” en la función plm (). El siguiente fragmento de código usa la muestra completa.

wage.within <- plm(lwage~exper+I(exper^2)+
                  tenure+I(tenure^2)+south+union,
                  data=nlspd, 
                  model="within")
tbl <- tidy(wage.within)
kable(tbl, digits=5, caption=
"Fixed effects using 'within' with full sample")

Fixed effects using ‘within’ with full sample
term	estimate	std.error	statistic	p.value
exper	0.04108	0.00662	6.20590	0.00000
I(exper^2)	-0.00041	0.00027	-1.49653	0.13463
tenure	0.01391	0.00328	4.24333	0.00002
I(tenure^2)	-0.00090	0.00021	-4.35357	0.00001
south	-0.01632	0.03615	-0.45153	0.65164
union	0.06370	0.01425	4.46879	0.00001

Ahora solo para una muestra de 10:

wage10.within <- plm(lwage~exper+I(exper^2)+
                  tenure+I(tenure^2)+union,
                  data=nls10, 
                  model="within")
tbl <- tidy(wage10.within)
kable(tbl, digits=5, caption=
  "Fixed effects using the 'within' model option for n=10")

Fixed effects using the ‘within’ model option for n=10
term	estimate	std.error	statistic	p.value
exper	0.23800	0.18776	1.26758	0.21332
I(exper^2)	-0.00819	0.00790	-1.03584	0.30738
tenure	-0.01235	0.03414	-0.36171	0.71974
I(tenure^2)	0.00230	0.00269	0.85407	0.39887
union	0.11354	0.15086	0.75263	0.45670

Los resutlados ahora obtenidos comparados con lso rpecedente meustran que los dos modelos son equivalentes (osea es igual a cuando utilizamos dummies)Mientras cuando comparamos con el modelo pooled donde podemos notar que inclueyendo los efectos fijos el efecto marginal de las varaibles disminuye notablemente. Podemos hacer un test para averiguar cual modelo se acomoda mejor por medio de la funcion pFtest():

kable(tidy(pFtest(wage.within, wage.pooled)), caption=
        "Fixed effects test: Ho:'No fixed effects'")

## Multiple parameters; naming those columns df1, df2

Fixed effects test: Ho:‘No fixed effects’
df1	df2	statistic	p.value	method	alternative
713	2858	15.18752	0	F test for individual effects	significant effects

La hipotesis nula de “no efectos fijos” es rechazada

Efectos aleatorios

El modelo de efectos aleatorios se basa en el modelo de efectos fijos reconociendo que, dado que los individuos en el panel se seleccionan al azar, sus características, medidas por la intersección $\beta_{1i}$β1i, también deben ser aleatorias. Por lo tanto, el modelo de efectos aleatorios asume la forma de la intersección como se indica en la Ecuación 3, donde $\bar \beta_1$ representa el promedio de la población y $u_i$ representa un término aleatorio específico de un individuo. Como en el caso de los efectos fijos, los efectos aleatorios también son invariantes en el tiempo.

\[ \beta_{1i} =\overline \beta_{1}+u_{i}...(3)\]

Sostituyendo en 2 obtenemos:

\[ y_{it}=\overline \beta_{1}+\beta_{2}x_{2it}+\nu_{it}...(4)\] El modelo de efectos aleatorio por lo tanto se distingue por la estructura de su termino de error. Los errores tienen media cero, verianza igual a $\sigma^2_u+ \sigma^2_e$, no correlacion entre los individuos y una covarianza a lo largo del tiempo igual al primer termino de varianza.

un aspecto improtante de este modelo es que la correlación temporal en los errores no disminuye con el tiempo.

\[\rho=corr(\nu_{it},\nu_{is})=\frac{\sigma_{u}^{2}}{\sigma_{u}^{2}+\sigma_{e}^{2}}...(6)\]

Probar efectos aleatorios equivale a probar la hipótesis de que no hay diferencias entre los individuos, lo que implica que la variable aleatoria específica del individuo tiene una varianza cero. Matematicamente:

\[H_{0}:\;\sigma_{u}^{2}=0,\;\;\;H_{A}:\sigma_{u}^{2}>0\]

La misma función que usamos para efectos fijos se puede usar para efectos aleatorios, pero estableciendo el argumento model = en ‘random’ y seleccionando el método random.method como una de cuatro posibilidades: “swar” (predeterminado), “amemiya”, ” walhus ”, o“ nerlove ”. La función de prueba de efectos aleatorios es plmtest (), que toma como argumento principal el modelo de agrupación.

wageReTest <- plmtest(wage.pooled, effect="individual")
kable(tidy(wageReTest), caption=
        "A random effects test for the wage equation")

A random effects test for the wage equation
statistic	p.value	method	alternative
62.12314	0	Lagrange Multiplier Test - (Honda) for balanced panels	significant effects

La tabla precedente muestra que la hipotesis nula de varianza cero en los errores specificos individuales es rechazada, por lo tanto podria haber heterogeneidad entre lso individuos significativa.

Los estimadores de efectos aleatorios son confiables bajo el supuesto de que las características individuales (heterogeneidad) son exógenas, es decir, son independientes con respecto a los regresores en la ecuación de efectos aleatorios. Podemos utilizar Hausman para probar eso con la hipótesis nula de que los efectos aleatorios individuales son exógenos. La función de prueba phtest () compara los modelos de efectos fijos y de efectos aleatorios; las siguientes líneas de código estiman el modelo de efectos aleatorios y realizan la prueba de endogeneidad de Hausman.

wage.random <- plm(lwage~educ+exper+I(exper^2)+
                  tenure+I(tenure^2)+black+south+union,
                  data=nlspd, random.method="swar",
                  model="random")
kable(tidy(wage.random), digits=4, caption=
      "The random effects results for the wage equation")

The random effects results for the wage equation
term	estimate	std.error	statistic	p.value
(Intercept)	0.5339	0.0799	6.6839	0.0000
educ	0.0733	0.0053	13.7417	0.0000
exper	0.0436	0.0064	6.8606	0.0000
I(exper^2)	-0.0006	0.0003	-2.1361	0.0327
tenure	0.0142	0.0032	4.4699	0.0000
I(tenure^2)	-0.0008	0.0002	-3.8790	0.0001
black	-0.1167	0.0302	-3.8643	0.0001
south	-0.0818	0.0224	-3.6505	0.0003
union	0.0802	0.0132	6.0724	0.0000

Ahora resumiendo:

kable(tidy(phtest(wage.within, wage.random)), caption=
 "Hausman endogeneity test for the random effects wage model")

Hausman endogeneity test for the random effects wage model
statistic	p.value	parameter	method	alternative
20.72521	0.0020552	6	Hausman Test	one model is inconsistent

El p value nos indica que se rechaza la hipótesis nula que dice que los efectos aleatorios individuales son exógenos, lo que hace que la ecuación de efectos aleatorios sea inconsistente. En este caso, el modelo de efectos fijos es la solución correcta. (El número de parámetros de la tabla 15.10 se proporciona solo para las variables que varían en el tiempo).

Sin embargo, el modelo de efectos fijos no admite variables invariantes en el tiempo como educ o black. Dado que el problema del modelo de efectos aleatorios es la endogeneidad, se pueden usar métodos de variables instrumentales cuando los regresores invariantes en el tiempo deben estar en el modelo.

El estimador de Hausman-Taylor utiliza variables instrumentales en un modelo de efectos aleatorios; asume cuatro categorías de regresores: exógenos variables en el tiempo, endógenos variables en el tiempo, exógenos invariantes en el tiempo y endógenos invariantes en el tiempo.

El número de variables que varían en el tiempo debe ser al menos igual al número de variables que no varían en el tiempo. En nuestro modelo salarial, suponga que exper, tenencia y unión son exógenos variables en el tiempo, sur es endógeno variable en el tiempo, negro es exógeno invariante en el tiempo y educ es endógeno invariante en el tiempo. La misma función plm () permite realizar la estimación de Hausman-Taylor estableciendo model = “ht”.

wage.HT <- plm(lwage~educ+exper+I(exper^2)+
      tenure+I(tenure^2)+black+south+union |
      exper+I(exper^2)+tenure+I(tenure^2)+union+black,
      data=nlspd, model="ht")

## Warning: uses of 'pht()' and 'plm(., model = "ht")' are discouraged,
## better use 'plm(., model = "random", random.method = "ht", inst.method =
## "baltagi"/"am"/"bms")' for Hausman-Taylor, Amemiya-MaCurdy, and Breusch-Mizon-
## Schmidt estimator

kable(tidy(wage.HT), digits=5, caption=
     "Hausman-Taylor estimates for the wage equation")

Hausman-Taylor estimates for the wage equation
term	estimate	std.error	statistic	p.value
(Intercept)	-0.75077	0.58624	-1.28066	0.20031
educ	0.17051	0.04446	3.83485	0.00013
exper	0.03991	0.00647	6.16382	0.00000
I(exper^2)	-0.00039	0.00027	-1.46222	0.14368
tenure	0.01433	0.00316	4.53388	0.00001
I(tenure^2)	-0.00085	0.00020	-4.31885	0.00002
black	-0.03591	0.06007	-0.59788	0.54992
south	-0.03171	0.03485	-0.91003	0.36281
union	0.07197	0.01345	5.34910	0.00000

Los resutlados son arrojados por la estimacion HT. Se nota que los cambios mas gradnes son determinados por educ y black.

Ejemplo: Función de costo de aerolíneas

Objetivo

Suponga que deseamos averiguar cómo se comporta el costo total (C) en relación con la producción (Q), el precio del combustible (PF) y el factor de carga (LF). En resumen, deseamos estimar la función de costos de la aerolínea.
La pregunta es: ¿cómo se toman en cuenta los efectos no observables, o heterogeneidad, para obtener estimaciones consistentes y eficientes de los parámetros de las variables de interés primordial, que son producción, precio del combustible y **factor de carga* en este caso?

Datos panel

datos<- read_excel("Table 16_1.xls", skip = 3)
head(datos,5)

## # A tibble: 5 x 10
##       I     T       C     Q     PF    LF  ln_C    ln_Q ln_PF  ln_LF
##   <dbl> <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl> <dbl>  <dbl>
## 1     1     1 1140640 0.953 106650 0.534  13.9 -0.0484  11.6 -0.626
## 2     1     2 1215690 0.987 110307 0.532  14.0 -0.0133  11.6 -0.630
## 3     1     3 1309570 1.09  110574 0.548  14.1  0.0880  11.6 -0.602
## 4     1     4 1511530 1.18  121974 0.541  14.2  0.162   11.7 -0.615
## 5     1     5 1676730 1.16  196606 0.591  14.3  0.149   12.2 -0.526

Los datos analizan los costos de seis líneas de aviación comercial de 1970 a 1984, para un total de 90 observaciones de datos de panel. Las variables se definen como sigue: $I$ identificación de la aerolínea; $T$ = identificación del año; $Q$ producción, como ingresos por milla por pasajero, un índice; $C$ costo total, en 1 000 dólares; $PF$ precio del combustible; y $LF$ factor de carga, la utilización promedio de la capacidad de la flotilla.

El modelo agrupado(The Pooled Model)

Es un modelo de la forma

\[ TC_{it}=\beta_{1}+\beta_2Q_{it}+\beta_3PF_{it}+\beta_4LF_{it}+u_{it} \] con $i=1,2,,,6$ y $t=1,2,3,...,15$

modelo.pooled <-plm(C~Q+PF+LF,model="pooling",data=datos)
kable(tidy(modelo.pooled), digits=3,caption="Regresión agrupada(Pooled model)")

Regresión agrupada(Pooled model)
term	estimate	std.error	statistic	p.value
(Intercept)	1158559.201	360592.720	3.213	0.002
Q	2026114.340	61806.945	32.781	0.000
PF	1.225	0.104	11.814	0.000
LF	-3065753.131	696327.318	-4.403	0.000

tbl<- tidy(coeftest(modelo.pooled, vcov=vcovHC(modelo.pooled, type="HC0",cluster="group")))
kable(tbl, digits=5, caption= "Modelo Pooled con errores estándar robustos cluster")

Modelo Pooled con errores estándar robustos cluster
term	estimate	std.error	statistic	p.value
(Intercept)	1158559.20109	408479.30313	2.83627	0.00569
Q	2026114.33991	79510.13402	25.48247	0.00000
PF	1.22535	0.34876	3.51345	0.00071
LF	-3065753.13068	1017050.92228	-3.01436	0.00338

Comentarios:

El problema principal de este modelo es que no distingue entre las diferentes aerolíneas ni indica si la respuesta de costo total a las variables explicativas a través del tiempo es la misma para todas las aerolíneas. En otras palabras, si agrupamos diferentes aerolíneas en diferentes periodos se oculta la heterogeneidad (individualidad o singularidad) que puede existir entre las aerolíneas.

Modelo de mínimos cuadrados con variable dicótoma (MCVD) de efectos fijos

Tenemos el siguiente modelo de (regresión) de efectos fijos (MEF):

\[\begin{equation} TC_{it}=\beta_{1i}+\beta_2Q_{it}+\beta_3PF_{it}+\beta_4LF_{it}+u_{it} \end{equation}\]

Utilizamos el subíndice $i$ en el término del intercepto para indicar que los interceptos de las seis aerolíneas pueden ser diferentes.
Las diferencias quizá se deban a características especiales de cada aerolínea, como el estilo de administración, la fi losofía de la empresa o el tipo de mercado que atiende cada aerolínea.
El término “efectos fijos” se debe a que, aunque el intercepto puede diferir entre los sujetos (en este caso las seis aerolíneas), el intercepto de cada entidad no varía con el tiempo, es decir, es invariante en el tiempo

#datosA12<-as.data.frame(datos[1:30,])
#datosA12$I<-as.factor(datosA12$I)
#datosA12$T<-as.factor(datosA12$T)

datos$I<-as.factor(datos$I)
datos$T<-as.factor(datos$T)

modeloA.1<-lm(C~Q,data = subset(datos,I=="1"))
modeloA.2<-lm(C~Q,data = subset(datos,I=="2"))
modeloA.3<-lm(C~Q,data = subset(datos,I=="3"))
modeloA.4<-lm(C~Q,data = subset(datos,I=="4"))
modeloA.5<-lm(C~Q,data = subset(datos,I=="5"))
modeloA.6<-lm(C~Q,data = subset(datos,I=="6"))

salida.modeloA.1<- data.frame(xtable(summary(modeloA.1))) #Rmodelo Aerolínea 1
kable(salida.modeloA.1, caption="Ecuación de Costo Aerolínea 1")

Ecuación de Costo Aerolínea 1
	Estimate	Std..Error	t.value	Pr…t..
(Intercept)	-2448181	572677.1	-4.274975	0.0009043
Q	3607341	396672.2	9.094010	0.0000005

salida.modeloA.2<- data.frame(xtable(summary(modeloA.2))) #modelo Aerolínea 2
kable(salida.modeloA.2, caption="Ecuación de Costo Aerolínea 2")

Ecuación de Costo Aerolínea 2
	Estimate	Std..Error	t.value	Pr…t..
(Intercept)	-1795170	465435.2	-3.85697	0.0019812
Q	3858635	439269.4	8.78421	0.0000008

salida.modeloA.3<- data.frame(xtable(summary(modeloA.3))) #modelo Aerolínea 3
kable(salida.modeloA.3, caption="Ecuación de Costo Aerolínea 3")

Ecuación de Costo Aerolínea 3
	Estimate	Std..Error	t.value	Pr…t..
(Intercept)	-557635.3	223226.7	-2.498067	0.0266859
Q	3105950.7	528651.7	5.875231	0.0000545

salida.modeloA.4<- data.frame(xtable(summary(modeloA.4))) #modelo Aerolínea 4
kable(salida.modeloA.4, caption="Ecuación de Costo Aerolínea 4")

Ecuación de Costo Aerolínea 4
	Estimate	Std..Error	t.value	Pr…t..
(Intercept)	-304044	86396.79	-3.519159	0.0037729
Q	4432153	374598.60	11.831739	0.0000000

salida.modeloA.5<- data.frame(xtable(summary(modeloA.5))) #modelo Aerolínea 4
kable(salida.modeloA.5, caption="Ecuación de Costo Aerolínea 5")

Ecuación de Costo Aerolínea 5
	Estimate	Std..Error	t.value	Pr…t..
(Intercept)	-106440.6	18925.89	-5.624074	0.0000828
Q	3511320.4	150811.65	23.282819	0.0000000

salida.modeloA.6<- data.frame(xtable(summary(modeloA.6))) #modelo Aerolínea 4
kable(salida.modeloA.6, caption="Ecuación de Costo Aerolínea 6")

Ecuación de Costo Aerolínea 6
	Estimate	Std..Error	t.value	Pr…t..
(Intercept)	-58733.76	7942.631	-7.394748	0.0000052
Q	3572168.78	60533.190	59.011739	0.0000000

#Calculamos interceptos específicos para las aerolíneas



interceptos <- c(coef(modeloA.1)["(Intercept)"],  #intercepto de Aerolínea 1
                coef(modeloA.2)["(Intercept)"],  #intercepto de Aerolínea 2
                coef(modeloA.3)["(Intercept)"], #intercepto de Aerolínea 3
                coef(modeloA.4)["(Intercept)"], #intercepto de Aerolínea 4
                coef(modeloA.5)["(Intercept)"], #intercepto de Aerolínea 5
                coef(modeloA.6)["(Intercept)"]) #intercepto de Aerolínea 6
              

rectas.df<- data.frame(interceptos = interceptos,
                       slopes =c(coef(modeloA.1)["Q"],coef(modeloA.2)["Q"]),
                      coef(modeloA.3)["Q"],coef(modeloA.4)["Q"],
                      coef(modeloA.5)["Q"],coef(modeloA.6)["Q"],
                  Aerolinea = levels(datos$I))

## Warning in data.frame(interceptos = interceptos, slopes = c(coef(modeloA.1)
## ["Q"], : row names were found from a short variable and have been discarded

qplot(x=Q, y=C,color= I,data = datos) + 
  geom_abline(aes(intercept = interceptos, 
                  slope = slopes, 
                  color = Aerolinea), data = rectas.df)+ylim(min(interceptos),max(datos$C))+
  xlim(min(datos$Q),max(datos$Q))

Gráfico: Diferencia entre el modelo de regresión agrupada y el modelo de efectos fijos

#Grafico de modelo individual de costos para dos aerolíneas  vs agrupado
#Calculamos interceptos específicos para ambas aerolíneas y el modelo agrupado

datosA12<-as.data.frame(datos[1:30,])
datosA12$I<-as.factor(datosA12$I)
datosA12$T<-as.factor(datosA12$T)


#Modelo agrupado para los datos de las aerolíneas 1 y 2 
modelo.pooled12<-plm(C~Q,model="pooling",data=datosA12)
kable(tidy(modelo.pooled12), digits=3,caption="Regresión agrupada(Pooled model)")

Regresión agrupada(Pooled model)
term	estimate	std.error	statistic	p.value
(Intercept)	-1269517	403181.7	-3.149	0.004
Q	3010361	318397.0	9.455	0.000

interceptos <- c(coef(modeloA.1)["(Intercept)"],  #intercepto de Aerolínea 1
                coef(modeloA.2)["(Intercept)"])  #intercepto de Aerolínea 2) 
              

rectas.df<- data.frame(interceptos = interceptos,
                       slopes =c(coef(modeloA.1)["Q"],coef(modeloA.2)["Q"]),
                  Aerolinea = levels(datos$I)[1:2])


qplot(x=Q, y=C,color= I,data = datosA12) + 
  geom_abline(aes(intercept = interceptos, 
                  slope = slopes, 
                  color = Aerolinea), data = rectas.df)+
 geom_abline(aes(intercept =coef(modelo.pooled12)["(Intercept)"], 
                  slope =coef(modelo.pooled12)["Q"], 
                  color = "POOLED"))+
  ylim(min(interceptos),max(datos$C))+xlim(min(datos$Q),max(datos$Q))

Observa cómo la regresión agrupada sesga la estimación de la pendiente al ignorar los efectos fijos

Ajuste de un Modelo de efectos fijos con variables dicotómicas

El modelo de efectos fijos con variables dicotómicas sería el siguiente:

\[C_{it}=\alpha_1+\alpha_2D_{2i}+\alpha_3D_{3i}+\alpha_4D_{4i}+\alpha_5D_{5i}+\alpha_6 D_{6i}+\beta_2Q_{it}+\beta_3PF_{it}+\beta_4LF_{it}+u_{it}\]

donde $D_{2i}=1$ si la observación corresponde a la aerolínea 2, y 0 en otro caso; $D_{3i} = 1$ si la observación es de la aerolínea 3, y 0 en otro caso; y así sucesivamente.

#Modelo de efectos fijos con variables dicotómicas
modelo.EF.dicotómicas<- lm(C~Q+PF+LF+factor(I),data=datos)
kable(tidy(modelo.EF.dicotómicas),digits=3,caption="Efectos fijos con variables dicotómicas")

Efectos fijos con variables dicotómicas
term	estimate	std.error	statistic	p.value
(Intercept)	-131235.989	350777.079	-0.374	0.709
Q	3319023.288	171354.071	19.369	0.000
PF	0.773	0.097	7.944	0.000
LF	-3797367.596	613773.136	-6.187	0.000
factor(I)2	601733.199	100895.702	5.964	0.000
factor(I)3	1337180.439	186171.009	7.183	0.000
factor(I)4	1777592.126	213162.873	8.339	0.000
factor(I)5	1828252.333	231229.679	7.907	0.000
factor(I)6	1706474.257	228300.874	7.475	0.000

Comentario:

Observamos que todos los coeficientes de los interceptos diferenciales son muy signifi cativos estadísticamente en lo individual, lo cual indica que tal vez las seis aerolíneas son heterogéneas y, por tanto, los resultados de la regresión agrupada pueden ser dudosos (sesgados).
Los valores de los coeficientes de las pendientes también son diferentes al modelo agrupado.
La prueba F restringida es una prueba formal para contratar los dos modelos: modelo agrupado vs efectos fijos vía variable dicotómica.

\[F = \frac{(SSR_{\text{restringido}} - SSR_{\text{completo}})/q}{SSR_{\text{completo}} / (n-k-1)}\] donde $q$ es el número de restricciones bajo la hipótesis nula y $k$ es el número de regresoras en el modelo no restringido(completo).

La hipótesis nula de la F restringida es que todos los interceptos diferenciales son iguales a cero.

linearHypothesis(modelo.EF.dicotómicas, c("factor(I)2=0","factor(I)3=0","factor(I)4=0",
                                          "factor(I)5=0","factor(I)6=0"))

## Linear hypothesis test
## 
## Hypothesis:
## factor(I)2 = 0
## factor(I)3 = 0
## factor(I)4 = 0
## factor(I)5 = 0
## factor(I)6 = 0
## 
## Model 1: restricted model
## Model 2: C ~ Q + PF + LF + factor(I)
## 
##   Res.Df           RSS Df     Sum of Sq      F          Pr(>F)    
## 1     86 6817712805793                                            
## 2     81 3586497349242  5 3231215456552 14.595 0.0000000003467 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comentario a la prueba F restringida:

En este ejemplo el valor F calculado de 5 gl para el numerador y 81 gl para el denominador es muy significativo estadísticamente. Por tanto, se rechaza la hipótesis nula de que todos los interceptos (diferenciales) son iguales a cero. Hay efectos fijos.
Si el valor F no fuera significativo estadísticamente, concluiríamos que no hay diferencias entre los interceptos de las seis aerolíneas y puedieras optar por un modelo agrupado.

Estimador de efectos fijos dentro del grupo (DG)

Efectos fijos vía OLS con valores corregidos por la media, o “sin media”.

# obtenemos los valores corregidos por la media (demeaned data)
datos_demeaned <- with(datos,
            data.frame(C_demeaned = C-ave(C,I),
                       Q_demeaned = Q-ave(Q,I),
                       PF_demeaned = PF- ave(PF,I),
                       LF_demeaned=LF-ave(LF, I)))
head(datos_demeaned)

##   C_demeaned Q_demeaned PF_demeaned  LF_demeaned
## 1 -1498363.3 -0.4574739   -355670.2 -0.062704667
## 2 -1423313.3 -0.4234739   -352013.2 -0.064863667
## 3 -1329433.3 -0.3182509   -351746.2 -0.049455667
## 4 -1127473.3 -0.2344509   -340346.2 -0.056345667
## 5  -962273.3 -0.2500609   -265714.2 -0.006024667
## 6  -815263.3 -0.2364709   -196711.2 -0.021774667

# Estimamos la regresión usando las variables corregidas por la media
modelo.EF.demeaned<-lm(C_demeaned~Q_demeaned+PF_demeaned+LF_demeaned-1,data=datos_demeaned)

kable(data.frame(xtable(summary(modelo.EF.demeaned))))

	Estimate	Std..Error	t.value
Q_demeaned	3319023.2883805	165339.7626030	20.073957
PF_demeaned	0.7730709	0.0939033	8.232630
LF_demeaned	-3797367.5961274	592230.4851160	-6.411976

Comentarios sobre la regresión

no tiene término de intercepto
variables que no cambian con el tiempo se eliminarían y no será posible capturar su efecto sobre la variable dependiente.

Ajuste del modelo de efectos aleatorios usando el paquete “plm”

modelo.within<- plm(C~Q+PF+LF,data=datos,model="within")
tbl<- tidy(modelo.within)
kable(tbl, digits=5, caption="Efectos fijos usando el estimador 'within'")

Efectos fijos usando el estimador ‘within’
term	estimate	std.error	statistic
Q	3319023.28838	171354.07093	19.36939
PF	0.77307	0.09732	7.94368
LF	-3797367.59613	613773.13572	-6.18692

# summary usando errores estándar robuestos
coeftest(modelo.within, vcov. = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##          Estimate     Std. Error t value             Pr(>|t|)    
## Q   3319023.28838   160943.38653 20.6223 < 0.0000000000000002 ***
## PF        0.77307        0.26989  2.8644              0.00532 ** 
## LF -3797367.59613  1643716.22257 -2.3102              0.02342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prueba de la F restrictiva usando pFtest del pkg “plm”

kable(tidy(pFtest(modelo.within, modelo.pooled)), caption=
        "Prueba  de efectos fijos: Ho:'No efectos fijos'")

## Multiple parameters; naming those columns df1, df2

Prueba de efectos fijos: Ho:‘No efectos fijos’
df1	df2	statistic	p.value	method	alternative
5	81	14.59521	0	F test for individual effects	significant effects

Comentario final: Observe que el intercepto estimado de cada aerolínea representa las características específicas de cada aerolínea, pero no podremos identificar estas características individualmente. Por consiguiente, el intercepto $\alpha_1$ de la aerolínea 1 representa la filosofía de la administración de esa aerolínea, la composición del consejo de administración, la personalidad del director general, el género del director general, etc. Todas estas características de heterogeneidad se integran al valor del intercepto.

Modelo de efectos aleatorios (MEFA)

La idea básica es comenzar con la ecuación

\[\begin{equation} TC_{it}=\beta_{1i}+\beta_2Q_{it}+\beta_3PF_{it}+\beta_4LF_{it}+u_{it} \end{equation}\]

En vez de considerar fija a $\beta_{1i}$, suponemos que es una variable aleatoria con un valor medio igual a $\beta_1$ (en este caso, sin subíndice $i$)
Además, el valor del intercepto para una empresa individual se expresa como:

\[\begin{equation} \beta_{1i}= \beta_1 + \varepsilon_i \end{equation}\]

donde $\varepsilon_i$ es un término de error aleatorio con valor medio igual a cero y varianza de $\sigma_{\varepsilon}^2$

-Lo que afirma el modelo de efectos aleatorios en esencia es que las seis empresas de la muestra se tomaron de un universo mucho más grande de este tipo de compañías, que tienen una media común para el intercepto $\beta_1$ y que las diferencias individuales en los valores del intercepto de cada empresa se reflejan en el término de error $\varepsilon_i$

Al sustituir las dos ecuaciones anteriores tenemos la siguiente expresión:

\[\begin{equation} TC_{it}=\beta_{1}+\beta_2Q_{it}+\beta_3PF_{it}+\beta_4LF_{it}+\varepsilon_i+u_{it} \end{equation}\]

\[\begin{equation} TC_{it}=\beta_{1}+\beta_2Q_{it}+\beta_3PF_{it}+\beta_4LF_{it}+w_i \end{equation}\]

con \[w_i= \varepsilon_i+u_{it}\]

El término de error compuesto $w_{it}$ consta de dos componentes:
- $\varepsilon_i$: componente de error de corte transversal o error específico del individuo,y
- $u_{it}$, la combinación del componente de error de series de tiempo y corte transversal, y que a veces se denomina término idiosincrásico porque varía en el corte transversal (es decir, el sujeto) así como en el tiempo.

Los supuestos comunes en los que se basa el Modelo de Efectos Aleatorios son:

$\varepsilon_i \sim \mathcal{N}(0,\sigma_{\varepsilon}^2)$
$u_{it}\sim\mathcal{N}(0,\sigma_{u}^2)$
$\mathbb{E}(\varepsilon_iu_{it})=0;\mathbb{E}(\varepsilon_i\varepsilon_j)=0,\quad \text{para} \quad i\neq j$
$\mathbb{E}(u_{it}u_{is})=0;\mathbb{E}(u_{it}u_{js})=0,\mathbb{E}(u_{it}u_{jt})=0\quad \text{para} \quad i\neq j; t\neq s$

Es muy importante observar que $w_{it}$ no está correlacionado con ninguna variable explicativa del modelo. Como $\varepsilon_i$ es un componente de $w_{it}$, es posible que el segundo esté correlacionado con las variables explicativas. Si en efecto es así, el modelo de efectos aleatorios producirá una estimación inconsistente de los coeficientes de regresión.

Como resultado de los supuestos establecidos se deriva que:

\[\mathbb{E}(w_{it})=0\] \[\mathbb{V}ar(w_{it})=\sigma_{\varepsilon}^{2}+\sigma_{u}^{2}\] Observa que si $\sigma_{\varepsilon}^2$ 0, no hay diferencia entre el modelo agrupado y el modelo de efectos aleatorios, en cuyo caso tan sólo se agrupan todas las observaciones (de corte transversal y de series de tiempo) y se lleva a cabo la regresión agrupada. Esto es válido porque en esta situación no hay efectos específicos del sujeto o porque todos se tomaron en cuenta en las variables explicativas.

El término de error $w_{it}$ es homoscedástico
Sin embargo, puede demostrarse que $w_{it}$ y $w_{is}$ ($t\neq s$) están correlacionados; es decir, los términos de error de una unidad de corte transversal dada en dos puntos en el tiempo están correlacionados. El coeficiente de correlación, $corr(w_{it}, w_{is})$, es el siguiente

\[\rho=corr(w_{it},w_{is})=\dfrac{\sigma_{\varepsilon}^2}{\sigma_{\varepsilon}^2+\sigma_u^2} \quad \text{con}\quad t\neq s\] Dos características especiales del coeficiente de correlación $\rho$.

Primera, para cualquier unidad de corte transversal dada, el valor de la correlación entre los términos de error en dos momentos sigue siendo el mismo, sin importar la distancia entre los dos periodos. L correlación no decrece con el tiempo.
Segunda, la estructura de correlación sigue siendo la misma para todas las unidades de corte transversal; es decir, es idéntica para todos los sujetos.
Se debe tomar en cuenta la estructura de correlación: El método más adecuado en este caso es el de mínimos cuadrados generalizados (MCG).

Prueba del multiplicador de Lagrange(Efectos aleatorios vs Modelo Pooled)

lagrangeTest<- plmtest(modelo.pooled, effect="individual",type = "bp")
lagrangeTest

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

kable(tidy(lagrangeTest), caption="Una prueba de efectos aleatorios vs modelo agrupado")

Una prueba de efectos aleatorios vs modelo agrupado
statistic	p.value	parameter	method	alternative
0.613087	0.4336279	1	Lagrange Multiplier Test - (Breusch-Pagan) for balanced panels	significant effects

Comentario:

Hipótesis nula: $\sigma_u^2=0$: No hay efectos individuales, equivale a un modelo pooled

Hipótesis alternativa:$\sigma_u^2\neq 0$: Hay efectos individuales, i.e, optar por modelo de efectos aleatorios.

La aplicación de la prueba BP produce un valor ji cuadrada de 0.61. El valor p de obtener un valor ji cuadrada de 0.61 o mayor es alrededor de 43%. Por consiguiente, no rechazamos la hipótesis nula. En otras palabras, el modelo de efectos aleatorios no es apropiado en el presente ejemplo.

modelo.efectos.aleatorios<- plm(C~Q+PF+LF,
                  data=datos, random.method="swar",
                  model="random")
kable(tidy(modelo.efectos.aleatorios), digits=4, caption=
      "Resultados del Modelo de Efectos Aleatorios para la función de Costo de las Aerolíneas")

Resultados del Modelo de Efectos Aleatorios para la función de Costo de las Aerolíneas
term	estimate	std.error	statistic	p.value
(Intercept)	1074292.9548	377467.9982	2.8461	0.0044
Q	2288587.9557	109493.7354	20.9015	0.0000
PF	1.1236	0.1034	10.8622	0.0000
LF	-3084994.0557	725679.8019	-4.2512	0.0000

Prueba de Hausman(Efectos aleatorios vs Efectos fijos)

La hipótesis nula en que se basa la prueba de Hausman es que los estimadores MEF y MCE no difieren considerablemente.
El estadístico de prueba desarrollado por Hausman tiene distribución $\chi^2$

\[\xi^H=(\hat{\beta}_{RE}-\hat{\beta}_{FE})\left[V(\hat{\beta}_{RE})-V(\hat{\beta}_{FE}\right]^{-1}(\hat{\beta}_{RE}-\hat{\beta}_{FE}) \sim \chi^2(K)\]

con $K=$ elementos comunes en $\beta$.

Si se rechaza la hipótesis nula, la conclusión es que el modelo de efectos aletorios no es apropiado porque es probable que los efectos aleatorios estén correlacionados con una o más regresoras. En este caso, se prefiere el modelo de efectos fijos.
Si no rechazas Hausman entonces te quedas con efectos aletorios.

#Prueba de Hausman para elegir entre Modelo de efectos aleatorios y Modelo de efectos fijos
phtest(modelo.within, modelo.efectos.aleatorios)

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

kable(tidy(phtest(modelo.within, modelo.efectos.aleatorios)), caption=
 "Prueba de Hausman")

Prueba de Hausman
statistic	p.value	parameter	method	alternative
60.86954	0	3	Hausman Test	one model is inconsistent

Comentario:

Hipótesis nula: los efectos fijos no están correlacionados con el error $\equiv$ modelo de efectos aleatorios consistente $\equiv$ optar por modelo de efectos aleatorios

Hipótesis aleternativa: los efectos fijos están correlacionados con el error $\equiv$ optar por modelo de efectos fijos

Como p-valor $<$ 0.05, rechazamos el Modelo de Efectos Aleatorios (MEFA) en favor del Modelo de Efectos Fijos(MEF). Esta conclusión se ve respaldada por la prueba del multiplicador de Lagrange de Breusch y Pagan pues igual se llegó a la conclusión de que el modelo de efectos aleatorios no es adecuado para el ejemplo de las aerolíneas.

Un cuadro resumen

Random Effects vs Fixed Effects

Referencias

Fuente de los datos

http://highered.mheducation.com/sites/0073375772/student_view0/data_sets.html

Laboratorio 9

Federico Daverio

Manejo de datos de Panel

Problema

Creaciòn de datos

Con el paquete reshape2

Ejemplo 1

El modelo pooled

Efectos fijos

Efectos aleatorios

Ejemplo: Función de costo de aerolíneas

Objetivo

Datos panel

El modelo agrupado(The Pooled Model)

Modelo de mínimos cuadrados con variable dicótoma (MCVD) de efectos fijos

Gráfico: Diferencia entre el modelo de regresión agrupada y el modelo de efectos fijos

Ajuste de un Modelo de efectos fijos con variables dicotómicas

Estimador de efectos fijos dentro del grupo (DG)

Efectos fijos vía OLS con valores corregidos por la media, o “sin media”.

Ajuste del modelo de efectos aleatorios usando el paquete “plm”

Prueba de la F restrictiva usando pFtest del pkg “plm”

Modelo de efectos aleatorios (MEFA)

Prueba del multiplicador de Lagrange(Efectos aleatorios vs Modelo Pooled)

Prueba de Hausman(Efectos aleatorios vs Efectos fijos)

Un cuadro resumen

Referencias

Fuente de los datos