#Para realizar el punto 3, lo primero que hacemos es crear los modelos lineales, cuadraticos y/o con interacción en relación a la cantidad de estimadores que especifica la función f(.)

Para ello lo primero que hacemos es importar los datos disponibles para posteriormente crear los modelos

d <- "C:/Users/MAEI/Documents/R"
setwd(d)

library(haven)
datos <- read_dta("datosfuncion.dta")
View(datos)

Antes de pasar a los modelos, veamos el comportamiento de las variables x1 y x2 en relación a la variable dependiente y

library(ggplot2)
x11()
ggplot (datos, aes (x = x1, y = y)) +
  geom_point (col = "blue", hight = 3)
## Warning: Ignoring unknown parameters: hight

Como podemo notar, la gráfica nos muestra que no hay alguna relación lineal directa, ni indirecta, ni no lineal para la variable x1

x11()
ggplot (datos, aes (x = x2, y = y)) +
  geom_point (col = "blue", hight = 3)
## Warning: Ignoring unknown parameters: hight

Como podemo notar, la gráfica nos muestra que no hay alguna relación lineal directa, ni indirecta, ni no lineal para la variable x2

Una vez obtenidos los datos, creamos los modelos:

Modelo Lineal

mod1 <- lm(y ~ x1 + x2, datos)
summary(mod1)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.633  -9.522  -2.088   9.694  18.522 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.1567     5.4212   3.718 0.000336 ***
## x1            2.9119     0.4145   7.025 2.96e-10 ***
## x2            1.3434     0.5597   2.400 0.018299 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.53 on 97 degrees of freedom
## Multiple R-squared:  0.3465, Adjusted R-squared:  0.333 
## F-statistic: 25.72 on 2 and 97 DF,  p-value: 1.095e-09

Modelo Lineal con interacción

mod2 <- lm(y ~ x1 + x2 + x1*x2, datos)
summary(mod2)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x1 * x2, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.835  -9.885  -1.922  10.469  19.137 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.0492    10.4803   3.726 0.000329 ***
## x1            1.1623     0.9298   1.250 0.214300    
## x2           -2.9683     2.1318  -1.392 0.167010    
## x1:x2         0.4072     0.1945   2.094 0.038938 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.34 on 96 degrees of freedom
## Multiple R-squared:  0.375,  Adjusted R-squared:  0.3555 
## F-statistic:  19.2 on 3 and 96 DF,  p-value: 7.794e-10

Modelo cuadratico

mod3 <- lm(y ~ x1 + x1^2 + x2 + x2^2, datos)
summary(mod3)
## 
## Call:
## lm(formula = y ~ x1 + x1^2 + x2 + x2^2, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.633  -9.522  -2.088   9.694  18.522 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.1567     5.4212   3.718 0.000336 ***
## x1            2.9119     0.4145   7.025 2.96e-10 ***
## x2            1.3434     0.5597   2.400 0.018299 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.53 on 97 degrees of freedom
## Multiple R-squared:  0.3465, Adjusted R-squared:  0.333 
## F-statistic: 25.72 on 2 and 97 DF,  p-value: 1.095e-09

Modelo cuadratico con interacción

mod4 <- lm(y ~ x1 + x1^2 + x2 + x2^2 + x2*x1, datos)
summary(mod4)
## 
## Call:
## lm(formula = y ~ x1 + x1^2 + x2 + x2^2 + x2 * x1, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.835  -9.885  -1.922  10.469  19.137 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.0492    10.4803   3.726 0.000329 ***
## x1            1.1623     0.9298   1.250 0.214300    
## x2           -2.9683     2.1318  -1.392 0.167010    
## x1:x2         0.4072     0.1945   2.094 0.038938 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.34 on 96 degrees of freedom
## Multiple R-squared:  0.375,  Adjusted R-squared:  0.3555 
## F-statistic:  19.2 on 3 and 96 DF,  p-value: 7.794e-10

Ahora analizamos hacemos el test de reset de Ramsey para saber cuál es el modelo mejor especificado

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
resettest(mod1, power = 2:3, type = "regressor")
## 
##  RESET test
## 
## data:  mod1
## RESET = 3.1838, df1 = 4, df2 = 93, p-value = 0.0169
resettest(mod2, power = 2:3, type = "regressor")
## 
##  RESET test
## 
## data:  mod2
## RESET = 2.4241, df1 = 4, df2 = 92, p-value = 0.05364
resettest(mod3, power = 2:3, type = "regressor")
## 
##  RESET test
## 
## data:  mod3
## RESET = 3.1838, df1 = 4, df2 = 93, p-value = 0.0169
resettest(mod4, power = 2:3, type = "regressor")
## 
##  RESET test
## 
## data:  mod4
## RESET = 2.4241, df1 = 4, df2 = 92, p-value = 0.05364

Con el test aplicado a los 4 modelos, podemos notar que con un p-valor mayor a 0.05 los unicos modelos que están aparentemente bien especificados, son el modelo lineal con interacción y el modelo cuadratico con interacción. El resto de modelos no pueden estar bien especificados por una omisión de una variable explicativa relevante, la forma funcional está incorrecta o errores de medición.

Notemos que el el modelo lineal con interacción y el modelo cuadratico tienen los mismos coeficientes al hacer la regresión:

library(car)
## Loading required package: carData
compareCoefs(mod2, mod4)
## Calls:
## 1: lm(formula = y ~ x1 + x2 + x1 * x2, data = datos)
## 2: lm(formula = y ~ x1 + x1^2 + x2 + x2^2 + x2 * x1, data = datos)
## 
##             Model 1 Model 2
## (Intercept)    39.0    39.0
## SE             10.5    10.5
##                            
## x1             1.16    1.16
## SE             0.93    0.93
##                            
## x2            -2.97   -2.97
## SE             2.13    2.13
##                            
## x1:x2         0.407   0.407
## SE            0.195   0.195
##

Debido a esto, para efectos practicos tomaremos el modelo lineal con interacción para analizar si hay cambio estructural en el modelo para saber si los betas se mantienen constantes a lo largo de toda la muestra

Cambio estructural

Para saber si hay cambio estructural aplicamos el test de Chow y lo hacemos manualmente

#Separando los datos en dos muestras 
datos.1 <- datos[1:50,]
datos.2 <- datos[51:100,]

mod1.1 <- lm(y~x1+x2+x1*x2, datos.1)
summary(mod1.1)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x1 * x2, data = datos.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7861 -1.8496  0.6207  2.5740  5.2611 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 32.68172    4.56004   7.167 5.12e-09 ***
## x1           0.96584    0.39109   2.470   0.0173 *  
## x2          -1.06832    0.89604  -1.192   0.2393    
## x1:x2        0.18147    0.08152   2.226   0.0309 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.409 on 46 degrees of freedom
## Multiple R-squared:  0.6781, Adjusted R-squared:  0.6571 
## F-statistic: 32.29 on 3 and 46 DF,  p-value: 2.183e-11
mod1.2 <- lm(y~x1+x2+x1*x2, datos.2)
summary(mod1.2)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x1 * x2, data = datos.2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.522 -2.847  1.548  3.419  6.096 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 23.35861    5.45036   4.286 9.21e-05 ***
## x1           3.80988    0.50233   7.584 1.22e-09 ***
## x2           0.66913    1.14682   0.583    0.562    
## x1:x2        0.02034    0.10546   0.193    0.848    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.157 on 46 degrees of freedom
## Multiple R-squared:  0.8879, Adjusted R-squared:  0.8806 
## F-statistic: 121.5 on 3 and 46 DF,  p-value: < 2.2e-16
#calculando la prueba de hipotesis para el cambio estructural
anova(mod2) #utilizo el anova para sacar la SCE del modelo
## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## x1         1  6076.9  6076.9 47.2645 6.255e-10 ***
## x2         1   766.5   766.5  5.9613   0.01645 *  
## x1:x2      1   563.5   563.5  4.3829   0.03894 *  
## Residuals 96 12342.9   128.6                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod1.1)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## x1         1 926.66  926.66 79.7164 1.317e-11 ***
## x2         1 141.94  141.94 12.2104  0.001063 ** 
## x1:x2      1  57.61   57.61  4.9556  0.030946 *  
## Residuals 46 534.73   11.62                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod1.2)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq  F value  Pr(>F)    
## x1         1 6128.3  6128.3 354.6706 < 2e-16 ***
## x2         1  168.2   168.2   9.7340 0.00312 ** 
## x1:x2      1    0.6     0.6   0.0372 0.84792    
## Residuals 46  794.8    17.3                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
SCE0 <- 12906.4
Fs <- ((12342.9-(534.73+794.8)/2)/((534.73+794.8)/(100-4)))
Fs
## [1] 843.231
Fp <- qf(0.05, 2, 96, lower.tail=FALSE)
Fp
## [1] 3.091191
#Como Fs > Fp se rechaza Ho y por lo tanto hay cambio estructural

Como lo podemos observa, el F teórico es menor al F del investigador por lo que se rechaza la hipotesis nula y se acepta la alterna, la cual dice que hay cambio estructural

Para verlo gráficamente utilicemos el test de cusum

library(strucchange)
## Loading required package: sandwich
cusum1<- efp(y~x1+x2+x1*x2,datos, type = "OLS-CUSUM")
x11()
plot(cusum1)

cusum3 <- efp(y~x1+x2+x1*x2,datos, type = "Rec-CUSUM")
plot(cusum3)

Como lo podemos notar en las dos gráficas, podemos reafirmar que hay cambio estructural en el modelo. Hay que tener en cuenta el cambio estructural en nuestro modelo ya que sabemos que debido a esto nuestros parametros poblacionales son diferentes a lo largo de la muestra, y por lo tanto la regresión del modelo puede no estar bien estimada.