length = 
      600;     
base  =  
      ReadList["C:\math\datos\datos1.dat", Table[Real, {4}], 
        length];
w = 
        Table[base[[i, 1]], {i, 1, length}];
x = 
        Table[base[[i, 2]], {i, 1, length}];
y = 
        Table[base[[i, 3]], {i, 1, length}];
z = Table[base[[i, 4]], {i, 1, length}];

"Variable W"
ListPlot[w, PlotJoined -> True]

    
"Variable X"
ListPlot[x, PlotJoined -> True]

    
"Variable Y"
ListPlot[y, PlotJoined -> True]

    
"Variable Z"
ListPlot[z, PlotJoined -> True]
"Se observa que tienen una distribucion normal con media y varianza \
constantes"

n = 
        Length[x];

fmv = L[σ_, (Subscript[β, 0]) _, (Subscript[β, 1]) _, (Subscript[β, 2]) _, (Subscript[β, 3]) _] = -(1/2)*Log[2 π] - 
            n*Log[σ] - (1/(2 σ)) (∑_(t = 1)^n)(((y[[t]] - (Subscript[β, 0])  - (Subscript[β, 1]) x[[t]] - (Subscript[β, 2]) z[[t]] - (Subscript[β, 3]) w[[t]]))^2);

dsigma = 
        Simplify[(∂_σ) fmv];

    db0 = Simplify[(∂_((Subscript[β, 0]))) fmv];

    db1 = Simplify[(∂_((Subscript[β, 1]))) fmv];

    db2 = Simplify[(∂_((Subscript[β, 2]))) fmv];

    db3 = Simplify[(∂_((Subscript[β, 3]))) fmv]

(*
Buscamos las raices,  
      que son en sì los estimadores de Máxima Verosimilitud*) 
raices = 
      FindRoot[{dsigma == 0, db0 == 0, db1 == 0, db2 == 0, 
          db3 == 
            0}, 
{σ,  .01}, 
{(Subscript[β, 0]), 3.4}, 
{(Subscript[β, 1]), 
          9.4}, 
{(Subscript[β, 2]), -80.4}, 
{(Subscript[β, 3]), 
          90.1}, 
MaxIterations -> 20000, 
        AccuracyGoal -> 10]

hess = {{(∂_σ) dsigma, (∂_((Subscript[β, 0]))) 
              dsigma, (∂_((Subscript[β, 1]))) 
              dsigma, (∂_((Subscript[β, 2]))) 
              dsigma, (∂_((Subscript[β, 3]))) 
              dsigma}, 
{(∂_σ) 
              db0, (∂_((Subscript[β, 0]))) 
              db0, (∂_((Subscript[β, 1]))) 
              db0, (∂_((Subscript[β, 2]))) 
              db0, (∂_((Subscript[β, 3]))) 
              db0}, 
{(∂_σ) 
              db1, (∂_((Subscript[β, 0]))) 
              db1, (∂_((Subscript[β, 1]))) 
              db1, (∂_((Subscript[β, 2]))) 
              db1, (∂_((Subscript[β, 3]))) 
              db1}, 
{(∂_σ) 
              db2, (∂_((Subscript[β, 0]))) 
              db2, (∂_((Subscript[β, 1]))) 
              db2, (∂_((Subscript[β, 2]))) 
              db2, (∂_((Subscript[β, 3]))) 
              db2}, 
{(∂_σ) 
              db3, (∂_((Subscript[β, 0]))) 
              db3, (∂_((Subscript[β, 1]))) 
              db3, (∂_((Subscript[β, 2]))) 
              db3, (∂_((Subscript[β, 3]))) db3}};

(*¿es una matriz simetrica?*) 
hess[[2, 1]] == hess[[1, 2]] /. 
      raices

    hess[[3, 1]] == hess[[1, 3]] /. raices

    hess[[3, 2]] == hess[[2, 3]] /. raices

    hess[[4, 1]] == hess[[1, 4]] /. raices

    hess[[4, 2]] == hess[[2, 4]] /. raices

    hess[[4, 3]] == hess[[3, 4]] /. raices

    hess[[5, 1]] == hess[[1, 5]] /. raices

    hess[[5, 2]] == hess[[2, 5]] /. raices

    hess[[5, 3]] == hess[[3, 5]] /. raices

    hess[[5, 4]] == hess[[4, 5]] /. raices

(*lo siguiente aplica si se trata de una matriz simétrica : 
1 Todos los valores propios son negativos  ->  
        definida negativa*) 
hessvalued = 
        hess /. raices;

    Eigenvalues[hessvalued]

    Det[hessvalued]

     (*si resulta negativo entonces se trata de un máximo*)

hess[[1, 
          1]] /. raices

    hess[[2, 2]] /. raices

    hess[[3, 3]] /. raices

    hess[[4, 4]] /. raices

    hess[[5, 5]] /. raices

hessvalued = 
        hess /. raices;

    Eigenvalues[hessvalued]