length = 
        600;    
base  =  
        ReadList["C:\math\datos\datos1.dat", Table[Real, {4}], 
          length];

    
x = 
        Table[base[[i, 2]], {i, 1, length}];
y = 
        Table[base[[i, 3]], {i, 1, length}];
"variable X"
ListPlot[x, PlotJoined -> True]
"variable Y"
ListPlot[y, PlotJoined -> True]

n = Length[y];

(*Función de verosimilitud en el caso Normal Bivariado = ((((1 - (ρ^2)))^(-1/2))/(2 π (Subscript[σ, 1]) (Subscript[σ, 2]))) 
          Exp[((-1)/(2 (1 - (ρ^2)))) ((((((x - (Subscript[μ, 1]))/(Subscript[σ, 1]))))^2) - 
                  2*ρ (((x - (Subscript[μ, 1]))/(Subscript[σ, 1]))) (((y - (Subscript[μ, 2]))/(Subscript[σ, 2]))) + (((((y - (Subscript[μ, 2]))/(Subscript[σ, 2]))))^2))]*)

L[ρ, (Subscript[μ, 2]) _, (Subscript[σ, 2]) _, (Subscript[μ, 1]) _, (Subscript[σ, 1]) _] = -(n/2)*Log[1 - (ρ^2)] - n*Log[2 π] - 
          n*Log[(Subscript[σ, 1])] - 
          n*Log[(Subscript[σ, 2])] - (1/(2 (1 - (ρ^2)))) (∑_(t = 1)^n)((((((x[[t]] - (Subscript[μ, 1]))/(Subscript[σ, 1]))))^2) - 
                  2*ρ (((x[[t]] - (Subscript[μ, 1]))/(Subscript[σ, 1]))) (((y[[t]] - (Subscript[μ, 2]))/(Subscript[σ, 2]))) + (((((y[[t]] - (Subscript[μ, 2]))/(Subscript[σ, 2]))))^2));

dmu1 = 
        Simplify[(∂_((Subscript[μ, 1]))) 
            L[ρ, (Subscript[μ, 2]) _, (Subscript[σ, 2]) _, (Subscript[μ, 1]) _, (Subscript[σ, 1]) _]];
dmu2 = 
        Simplify[(∂_((Subscript[μ, 2]))) 
            L[ρ, (Subscript[μ, 2]) _, (Subscript[σ, 2]) _, (Subscript[μ, 1]) _, (Subscript[σ, 1]) _]];
dsigma1 = 
        Simplify[(∂_((Subscript[σ, 1]))) 
            L[ρ, (Subscript[μ, 2]) _, (Subscript[σ, 2]) _, (Subscript[μ, 1]) _, (Subscript[σ, 1]) _]];
dsigma2 = 
        Simplify[(∂_((Subscript[σ, 2]))) 
            L[ρ, (Subscript[μ, 2]) _, (Subscript[σ, 2]) _, (Subscript[μ, 1]) _, (Subscript[σ, 1]) _]];
drho = 
        Simplify[(∂_ρ) 
            L[ρ, (Subscript[μ, 2]) _, (Subscript[σ, 2]) _, (Subscript[μ, 1]) _, (Subscript[σ, 1]) _]];

FindRoot[{dmu1 == 0, dmu2 == 0, dsigma1 == 0, dsigma2 == 0, 
        drho == 0}, 
{(Subscript[μ, 1]), 
        1}, 
{(Subscript[μ, 2]), 
        0}, 
{(Subscript[σ, 1]), 
        1}, 
{(Subscript[σ, 2]), 
        1}, 
{ρ, -0.14}, 
MaxIterations -> 20000, 
      AccuracyGoal -> 30]

Clear[x]
(Subscript[μ, 2]) + ((ρ  (Subscript[σ, 2]))/(Subscript[σ, 1])) 
          X - ((ρ (Subscript[σ, 2]) (Subscript[μ, 1]))/(Subscript[σ, 1])) /. {%%}