I estimated linear regression by using maximum likelihood in R. I did regression, put the data into matrices for the MLE procedure, and estimated the model. I used the following code and I got the following error. I guess that I need to fix the following: (c(15,0.9,13,15,0.9,13,15,0.9) that sets the starting values for β0, β1, β2, β3, β4, β5, β6,σ2. In addition, the code for t-statistics has an error. What should I do to fix it?
> x1 <- cbind(1,as.matrix(data$V2))
> x2 <- cbind(1,as.matrix(data$V3))
> x3 <- cbind(1,as.matrix(data$V4))
> x4 <- cbind(1,as.matrix(data$V5))
> x5 <- cbind(1,as.matrix(data$V6))
> x6 <- cbind(1,as.matrix(data$V7))
> y <- as.matrix(data$V1)
> ones <- x1[,1]
> ones2 <- x2[,1]
> ones3 <- x3[,1]
> ones4 <- x4[,1]
> ones5 <- x5[,1]
> ones6 <- x6[,1]
> K <- ncol(x1)
> K
[1] 2
> K1 <- K + 1
> Ka <- ncol(x2)
> Ka
[1] 2
> K2 <- Ka + 1
> Kb <- ncol(x3)
> Kb
[1] 2
> K3 <- Kb + 1
> Kc <- ncol(x4)
> Kc
[1] 2
> K4 <- Kc + 1
> Kd <- ncol(x5)
> Kd
[1] 2
> K5 <- Kd + 1
> Ke <- ncol(x6)
> Ke
[1] 2
> K6 <- Ke + 1
> n1 <- nrow(x1)
> n1
[1] 286
> n2 <- nrow(x2)
> n2
[1] 286
> n3 <- nrow(x3)
> n3
[1] 286
> n4 <- nrow(x4)
> n4
[1] 286
> n5 <- nrow(x5)
> n5
[1] 286
> n6 <- nrow(x6)
> n6
[1] 286
> llik.regress <- function(par,X1,X2, X3, X4, X5, X6, Y) {
+
+ Y <- as.vector(y)
+
+ X1 <- as.matrix(x1)
+
+ X2 <- as.matrix(x2)
+
+ X3 <- as.matrix(x3)
+
+ X4 <- as.matrix(x4)
+
+ X5 <- as.matrix(x5)
+
+ X6 <- as.matrix(x6)
+
+ xbeta1 <- X1%*%par[1:K]
+
+ xbeta2 <- X2%*%par[1:K]
+
+ xbeta3 <- X3%*%par[1:K]
+
+ xbeta4 <- X4%*%par[1:K]
+
+ xbeta5 <- X5%*%par[1:K]
+
+ xbeta6 <- X6%*%par[1:K]
+
+ Sig <- par[K1:K1]
+
+ sum(-(1/2)*log(2*pi)-(1/2)*log(Sig^2)-(1/(2*Sig^2))*((y-xbeta1)^2)*((y-xbeta2)^2)*((y-xbeta3)^2)*((y-xbeta4)^2)*((y-xbeta5)^2)*((y-xbeta6)^2))
+
+ }
> llik.regress
function(par,X1,X2, X3, X4, X5, X6, Y) {
Y <- as.vector(y)
X1 <- as.matrix(x1)
X2 <- as.matrix(x2)
X3 <- as.matrix(x3)
X4 <- as.matrix(x4)
X5 <- as.matrix(x5)
X6 <- as.matrix(x6)
xbeta1 <- X1%*%par[1:K]
xbeta2 <- X2%*%par[1:K]
xbeta3 <- X3%*%par[1:K]
xbeta4 <- X4%*%par[1:K]
xbeta5 <- X5%*%par[1:K]
xbeta6 <- X6%*%par[1:K]
Sig <- par[K1:K1]
sum(-(1/2)*log(2*pi)-(1/2)*log(Sig^2)-(1/(2*Sig^2))*((y-xbeta1)^2)*((y-xbeta2)^2)*((y-xbeta3)^2)*((y-xbeta4)^2)*((y-xbeta5)^2)*((y-xbeta6)^2))
}
> model <- optim(c(15,0.9,13,15,0.9,13,15,0.9),llik.regress, method = "BFGS", control = list(trace=6,maxit=100,fnscale = -1),
+ hessian = TRUE)
initial value 5394100715040953038084886.000000
iter 10 value 9717560928.044348
iter 20 value 35767262.281884
iter 30 value 143650.399866
final value 26963.930000
converged
> model
$par
[1] 3.967752e+00 2.259643e-04 1.815630e+03 1.500000e+01 9.000000e-01
[6] 1.300000e+01 1.500000e+01 9.000000e-01
$value
[1] -26963.93
$counts
function gradient
371 35
$convergence
[1] 0
$message
NULL
$hessian
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] -1.073836e+04 -5.439814e+08 7.809896e+00 0 0 0 0 0
[2,] -5.439814e+08 -2.792297e+13 -1.332402e+05 0 0 0 0 0
[3,] 7.809896e+00 -1.332402e+05 -4.460435e-02 0 0 0 0 0
[4,] 0.000000e+00 0.000000e+00 0.000000e+00 0 0 0 0 0
[5,] 0.000000e+00 0.000000e+00 0.000000e+00 0 0 0 0 0
[6,] 0.000000e+00 0.000000e+00 0.000000e+00 0 0 0 0 0
[7,] 0.000000e+00 0.000000e+00 0.000000e+00 0 0 0 0 0
[8,] 0.000000e+00 0.000000e+00 0.000000e+00 0 0 0 0 0
> v <- -solve(model$hessian)
Error in solve.default(model$hessian) :
Lapack routine dgesv: system is exactly singular: U[4,4] = 0
> v
Error: object 'v' not found
> se <- sqrt( diag(v))
Error in diag(v) : object 'v' not found
>
> se
Error: object 'se' not found
>
> t<-model$par/se
Error: object 'se' not found
> modelval1<-2*(1-pt(abs(t),nrow(X1)-ncol(X1)))
Error in nrow(X1) : object 'X1' not found
> modelval2<-2*(1-pt(abs(t),nrow(X2)-ncol(X2)))
Error in nrow(X2) : object 'X2' not found
> modelval3<-2*(1-pt(abs(t),nrow(X3)-ncol(X3)))
Error in nrow(X3) : object 'X3' not found
> modelval4<-2*(1-pt(abs(t),nrow(X4)-ncol(X4)))
Error in nrow(X4) : object 'X4' not found
> modelval5<-2*(1-pt(abs(t),nrow(X5)-ncol(X5)))
Error in nrow(X5) : object 'X5' not found
> modelval6<-2*(1-pt(abs(t),nrow(X6)-ncol(X6)))
Error in nrow(X6) : object 'X6' not found
> results<-cbind(model$par,se,t,modelval1,modelval2,modelval3,modelval4,modelval5,modelval6)
Error in cbind(model$par, se, t, modelval1, modelval2, modelval3, modelval4, :
object 'se' not found
> results(colnames)<-c("b","se","t","p")
Error in results(colnames) <- c("b", "se", "t", "p") :
could not find function "results<-"
> results(rownames)<-c("Const","X1","X2","X3","X4","X5","X6","Sigma2")
Error in results(rownames) <- c("Const", "X1", "X2", "X3", "X4", "X5", :
could not find function "results<-"
> print(results,digits=3)
Error in print(results, digits = 3) : object 'results' not found