Programe esa iteración en una rutina en R en la que entre la matriz \(\Sigma\), pruebe que esta matriz es definida positiva y retorne un codigo de error que indique la razón por la cual no se pudo realizar la descomposición y cero si se pudo realizar, en cuyo caso retorna L
choleskinsy<-function(A){
if(matrixNormal::is.square.matrix(A)==FALSE){stop("la matriz no es cuadrada")}
if(matrixNormal::is.symmetric.matrix(A)==FALSE){stop("la matriz no es simetrica")}
if(matrixNormal::is.positive.definite(A)==FALSE){stop("la matriz no es definida positiva")}
L<-matrix(0,nrow=nrow(A),ncol=ncol(A)) #creamos matriz de ceros para llenar
for(i in 1:nrow(A)){ #para cada fila
for(j in 1:i){ #resolvemos por columna
s<-0 #el valor de s comienza en cero para la primera columna según lo resuelto en el literal a)
if(j>1){ #Si el cálculo es en columnas despues de la primera
for(k in 1:(j-1)){ #establecemos el rango de la sumatoria desde k hasta j-1
s<-s+L[i,k]*L[j,k] #de esta forma el loop actualiza el valor de s añadiendo el término de cada columna adicional a la sumatoria
}}
if(i==j){#calculamos las diagonales
L[i,j]<-sqrt(A[i,j]-s)
}else {#calculamos los demás elementos de la matriz
L[i,j]<-(1/L[j,j])*(A[i,j]-s)
}
}
}
return(L)
} Encuentre los valores de una matriz simétrica definida positiva y léala o generela en R para verificar que su rutina hace bien el trabajo Comprobamos que funcione generando una matriz definida positiva
library(clusterGeneration)## Warning: package 'clusterGeneration' was built under R version 3.6.3## Loading required package: MASS## Warning: package 'MASS' was built under R version 3.6.3A<-genPositiveDefMat("unifcorrmat", dim=5)Empleamos la descomposición que creamos
L1<-choleskinsy(A$Sigma)
L1## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.1406478 0.0000000 0.000000 0.000000 0.000000
## [2,] -0.7233690 1.4387453 0.000000 0.000000 0.000000
## [3,] 0.8377504 0.5629948 2.096888 0.000000 0.000000
## [4,] -0.6745345 -0.1375418 2.180437 1.315164 0.000000
## [5,] 0.4760146 -1.1653550 -1.291946 1.523807 1.978249Ahora comparamos el resultado empleando la funcion chol de R. La función chol genera por defecto una matriz triangular superior, por lo que la trasponemos para que sean comparables los resultados
L2<-chol(A$Sigma)
L2<-t(L2)
L2## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.1406478 0.0000000 0.000000 0.000000 0.000000
## [2,] -0.7233690 1.4387453 0.000000 0.000000 0.000000
## [3,] 0.8377504 0.5629948 2.096888 0.000000 0.000000
## [4,] -0.6745345 -0.1375418 2.180437 1.315164 0.000000
## [5,] 0.4760146 -1.1653550 -1.291946 1.523807 1.978249Los resultados son los mismos!
Realice una simulación (todos los parámetros \(\sigma^2\), \(\delta_0\), \(\delta_1\), \(\beta_{2×1}\), \(\rho\), T = 1000) conocidos
# En primer lugar creamos los insumos del modelo
s_e2<-0.6
d0<-0.7
d1<-0.5
b_xx<-c(0.5,0.8)
rho<-0.6
m_xx<-c(15,12)
s_xx<-matrix(c(3,0,0,5),ncol=2)
# Creamos la simulación de los datos
simulacion_e2<-function(s,d0,d1,b,r,T){
X<-MASS::mvrnorm(T,m_xx,s_xx)
Z<-rpois(T,12)
p<-diag(1,nrow=T,T)
matriz_rho<-function(a,b){
for(i in 2:nrow(a)){
a[i, i-1]<-b
a[i-1, i]<-b
}
return(a)
}
Sig1<-matriz_rho(p,r)
Sig2<-(d0+d1*Z)%*%Sig1
Sig3<-matrix(0,nrow=T,ncol=T)
diag(Sig3)<-t(Sig2)
Sig4<-s*Sig3
e<-MASS::mvrnorm(n=1,mu=rep(0,T),Sig4)
Y<-X%*%b+e
datos_e2<-list()
datos_e2$X<-X
datos_e2$Y<-Y
datos_e2$e<-e
datos_e2$Sig<-Sig3
return(datos_e2)
}
# obtenemos los datos de la simulación
dat2<-simulacion_e2(s_e2,d0,d1,b_xx,rho,1000)Use su nueva rutina para hallar el modelo estrella \((Y^*,X^*)\), bajo el supuesto que conoce los parametros \(\sigma^2\) y \(\Sigma\) y halle el estimador OLS de \(\beta\) con base en los datos del modelo estrella
Mod_estrella<-function(dat2){
L<-choleskinsy(dat2$Sig)
C<-solve(L)
Y2<-C%*%dat2$Y
X2<-C%*%dat2$X
dat3<-data.frame(Y2,X2)
mod_e<-lm(Y2~X1+X2+0,data=dat3)
return(mod_e)
}
m2_4<-Mod_estrella(dat2)
summary(m2_4)##
## Call:
## lm(formula = Y2 ~ X1 + X2 + 0, data = dat3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.25031 -0.55714 -0.00856 0.54215 2.28640
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## X1 0.49933 0.02900 17.22 <2e-16 ***
## X2 0.78224 0.03595 21.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7804 on 998 degrees of freedom
## Multiple R-squared: 0.9709, Adjusted R-squared: 0.9709
## F-statistic: 1.667e+04 on 2 and 998 DF, p-value: < 2.2e-16Halle el estimador OLS de \(\beta\) con los datos originales y compare con los resultados del paso anterior. Que pasa si ahora simula T = 5000 o T = 10000
dat4<-data.frame(dat2$Y,dat2$X)
m2_5<-lm(dat2.Y~X1+X2+0,dat4)
summary(m2_5)##
## Call:
## lm(formula = dat2.Y ~ X1 + X2 + 0, data = dat4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.4084 -2.0648 -0.0241 2.0459 9.2491
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## X1 0.50468 0.02951 17.1 <2e-16 ***
## X2 0.77547 0.03657 21.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.998 on 998 degrees of freedom
## Multiple R-squared: 0.9699, Adjusted R-squared: 0.9698
## F-statistic: 1.605e+04 on 2 and 998 DF, p-value: < 2.2e-16La simulación con 10000 datos puede tomar mucho tiempo, por lo que solamente se presenta el resultado con 5000. En todo caso, el estimador se aproxima más al real a medida que incrementa el tamaño de muestra
dat5<-simulacion_e2(s_e2,d0,d1,b_xx,rho,5000)
m2_6<-Mod_estrella(dat5)
summary(m2_6)##
## Call:
## lm(formula = Y2 ~ X1 + X2 + 0, data = dat3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.74027 -0.51252 0.00424 0.51701 2.59234
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## X1 0.51666 0.01271 40.63 <2e-16 ***
## X2 0.77797 0.01579 49.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7669 on 4998 degrees of freedom
## Multiple R-squared: 0.9721, Adjusted R-squared: 0.9721
## F-statistic: 8.715e+04 on 2 and 4998 DF, p-value: < 2.2e-16