R Under development (unstable) (2023-08-27 r85021 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> library("mlt")
Loading required package: basefun
Loading required package: variables
> library("truncreg")
Loading required package: maxLik
Loading required package: miscTools

Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.

If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
> ### make M1mac happy
> options(digits = 4)
> 
> 
> ## Left-truncated 
> set.seed(29)
> n <- 1000	
> x <- runif(n, max = 2 * pi)
> y <- rnorm(n, 2*x + 1, .25)
> d <- data.frame(y = y, x = x)
> ## truncated response
> tr <- 2.5
> d$yt <- ifelse(d$y > tr, d$y, NA)
> d$trunc_left <- tr
> 
> tmod <- truncreg(yt ~ x, data = d, point = tr, direction = "left")
> coef(tmod)
(Intercept)           x       sigma 
     1.0193      1.9918      0.2575 
> logLik(tmod)
'log Lik.' -33.99 (df=3)
> vcov(tmod)
            (Intercept)          x      sigma
(Intercept)   4.972e-04 -1.147e-04 -7.445e-06
x            -1.147e-04  3.146e-05  1.637e-06
sigma        -7.445e-06  1.637e-06  3.941e-05
> 
> ## MLT
> fm <- as.formula("~ x")
> yb <- polynomial_basis(numeric_var("yt", support = range(d$yt, na.rm = TRUE)), 
+                        coef = c(TRUE, TRUE), ci = c(-Inf, 0))
> m <- ctm(yb, shifting = fm, todistr = "Normal", data = d)
> dfull <- d[!is.na(d$yt),,drop = FALSE]
> dfull$yt <- R(dfull$yt, tleft = dfull$trunc_left)
> mltmod <- mlt(m, data = dfull)
> (cf <- coef(mltmod))
(Intercept)          yt           x 
     -3.958       3.883      -7.734 
> c(-cf[1] / cf[2], -cf[3] / cf[2], 1 / cf[2])
(Intercept)           x          yt 
     1.0193      1.9918      0.2575 
> logLik(mltmod)
'log Lik.' -33.99 (df=3)
> vcov(mltmod)
            (Intercept)        yt        x
(Intercept)    0.017689 -0.009565  0.01723
yt            -0.009565  0.008957 -0.01774
x              0.017226 -0.017744  0.03563
> 
> #library("numDeriv")
> 
> solve(numDeriv::hessian(mltmod$loglik, coef(mltmod), 
+                         weights = weights(mltmod)))
          [,1]      [,2]     [,3]
[1,]  0.017689 -0.009565  0.01723
[2,] -0.009565  0.008957 -0.01774
[3,]  0.017226 -0.017744  0.03563
> 
> 
> ## right-truncated
> 
> set.seed(29)
> n <- 1000	
> x <- runif(n, max = 2 * pi)
> y <- rnorm(n, 2*x + 1, .25)
> d <- data.frame(y = y, x = x)
> ## truncated response
> tr <- 11
> d$yt <- ifelse(d$y < tr, d$y, NA)
> d$trunc_right <- tr
> 
> tmod <- truncreg(yt ~ x, data = d, point = tr, direction = "right")
> coef(tmod)
(Intercept)           x       sigma 
      1.000       1.995       0.264 
> logLik(tmod)
'log Lik.' -50.52 (df=3)
> vcov(tmod)
            (Intercept)          x      sigma
(Intercept)   3.312e-04 -1.022e-04 -2.841e-06
x            -1.022e-04  4.353e-05  1.855e-06
sigma        -2.841e-06  1.855e-06  4.526e-05
> 
> ## MLT
> fm <- as.formula("~ x")
> yb <- polynomial_basis(numeric_var("yt", support = range(d$yt, na.rm = TRUE)), 
+                        coef = c(TRUE, TRUE), ci = c(-Inf, 0))
> m <- ctm(yb, shifting = fm, todistr = "Normal", data = d)
> dfull <- d[!is.na(d$yt),,drop = FALSE]
> dfull$yt <- R(dfull$yt, tright = dfull$trunc_right)
> mltmod <- mlt(m, data = dfull)
> (cf <- coef(mltmod))
(Intercept)          yt           x 
     -3.789       3.788      -7.557 
> c(-cf[1] / cf[2], -cf[3] / cf[2], 1 / cf[2])
(Intercept)           x          yt 
      1.000       1.995       0.264 
> logLik(mltmod)
'log Lik.' -50.52 (df=3)
> vcov(mltmod)
            (Intercept)        yt        x
(Intercept)    0.014383 -0.009473  0.01733
yt            -0.009473  0.009316 -0.01849
x              0.017333 -0.018485  0.03730
> 
> #library("numDeriv")
> 
> solve(numDeriv::hessian(mltmod$loglik, coef(mltmod), 
+                         weights = weights(mltmod)))
          [,1]      [,2]     [,3]
[1,]  0.014245 -0.009332  0.01705
[2,] -0.009332  0.009171 -0.01820
[3,]  0.017050 -0.018196  0.03673
> 
> 
> 
> proc.time()
   user  system elapsed 
   2.07    0.42    2.45