# linear regression by using maximum likelihood in R.

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
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

> se <- sqrt( diag(v))
>
> se
>

> t<-model\$par/se
> modelval1<-2*(1-pt(abs(t),nrow(X1)-ncol(X1)))
> modelval2<-2*(1-pt(abs(t),nrow(X2)-ncol(X2)))
> modelval3<-2*(1-pt(abs(t),nrow(X3)-ncol(X3)))
> modelval4<-2*(1-pt(abs(t),nrow(X4)-ncol(X4)))
> modelval5<-2*(1-pt(abs(t),nrow(X5)-ncol(X5)))
> modelval6<-2*(1-pt(abs(t),nrow(X6)-ncol(X6)))