Clear[
      base]
  
Clear[lenght]
lenght = 237;
base  =  
        ReadList["C:\math\datos\pttopeoc.dat", Table[Real, {2}], 
          lenght];

pt = 
        Table[base[[i, 1]], {i, 1, lenght}];
topeoc = Table[base[[i, 2]], {i, 1, lenght}];

"Productividad laboral"
ListPlot[pt, PlotJoined -> False]
"Tope de ocupacion"
ListPlot[topeoc, PlotJoined -> False]
"Diagrama de dispersión"
ListPlot[ base, 
      PlotStyle -> {PointSize[0.02], Hue[0.6]}]

(*MODELO EXPONENCIAL GUMBEL : 
gfxy = ((1 + θ x) (1 + θ y) - θ) 
          Exp[-(x + 
                  y + θ x y)]
*) 
n = 3

    fmvexasprod = (∏_(i = 1)^n)(((1 + θ topeoc[[i]]) (1 + θ pt[[i]]) - θ) 
            Exp[-(topeoc[[i]] + 
                    pt[[i]] + θ*
                      topeoc[[i]] pt[[i]])])

fmvex = L[θ_] = (∑_(t = 1)^n)(Log[
                1 + θ topeoc[[t]] + θ pt[[t]] + (θ^2) topeoc[[t]] pt[[t]]  - θ] - topeoc[[t]] - 
              pt[[t]] - 
              topeoc[[t]] pt[[t]] θ)

dtheta = 
        Simplify[(∂_θ) fmvex];
raices = 
      FindRoot[{dtheta == 0}, 
{θ, 2.6,  .00001, 
          100}, 
MaxIterations -> 20000, 
        AccuracyGoal -> 10]
ddtheta = 
        Simplify[(∂_θ) dtheta];
ddtheta /. {raices}

Plot[fmvexasprod, {θ, 0,  .00001}]

fmvexasprod /. θ -> -0.01
fmvexasprod /. θ -> 0
fmvexasprod /. θ -> 0.01
fmvexasprod /. θ -> 0.05
fmvexasprod /. θ -> 0.5
fmvexasprod /. θ -> 1

?FindMaximum

li = 40;
ls = 80;

buscar = 10;

FindMaximum[
      fmvexasprod, {θ, 
        buscar}]
FindMaximum[
      fmvexasprod, {θ, buscar}, 
      Method -> 
        ConjugateGradient]
FindMaximum[
      fmvexasprod, {θ, buscar}, 
      Method -> 
        LevenbergMarquardt]
FindMaximum[
      fmvexasprod, {θ, buscar}, 
      Method -> Newton]
FindMaximum[fmvexasprod, {θ, buscar}, 
        Method -> QuasiNewton]

(*Metodos numericos:"NelderMead","DifferentialEvolution","SimulatedAnnealing"and "RandomSearch" *)
NMaximize[
      fmvexasprod, {θ, li, ls}, 
      Method -> NelderMead]
NMaximize[
      fmvexasprod, {θ, li, ls}, 
      Method -> 
        DifferentialEvolution]
NMaximize[
      fmvexasprod, {θ, li, ls}, 
      Method -> 
        SimulatedAnnealing]
NMaximize[
      fmvexasprod, {θ, li, ls}, 
      Method -> RandomSearch]

fmvexasprod /. θ -> -3778.036017112422`

(*ConjugateGradient","Gradient","LevenbergMarquardt","Newton"and"QuasiNewton"* *)

Limit[
      fmvexasprod, {θ -> Infinity}]
Limit[fmvexasprod, {θ -> 0}]
Limit[fmvexasprod, {θ -> 100}]
Limit[
      fmvexasprod, {θ -> -Infinity}]

Plot[fmvexasprod, {θ, 0, 0.00000000005}]

fmvexasprod /. θ -> 0

Min[pt]
Min[topeoc]