R Under development (unstable) (2026-03-22 r89674 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2026 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( "micEconDistRay" )

Please cite the 'micEconDistRay' package in publications as:
Price, J.J. and Henningsen, A. (2023): A Ray-Based Input Distance Function to Model Zero-Valued Output Quantities: Derivation and an Empirical Application. Journal of Productivity Analysis 60, p. 179-188. DOI: 10.1007/s11123-023-00684-1. URL: https://doi.org/10.1007/s11123-023-00684-1.

If you have questions, suggestions, or comments regarding the 'micEconDistRay' package, please use 'Issues' at the package's GitHub site:
https://github.com/micEcon/micEconDistRay/issues
> library( "quadprog" )
> 
> # load and prepare data set 
> data( appleProdFr86, package = "micEcon" )
> appleProdFr86$qCap <- appleProdFr86$vCap / appleProdFr86$pCap
> appleProdFr86$qLab <- appleProdFr86$vLab / appleProdFr86$pLab
> appleProdFr86$qMat <- appleProdFr86$vMat / appleProdFr86$pMat
> 
> # inputs
> xNames <- c( "qCap", "qLab", "qMat" )
> # outputs 
> yNames <- c( "qApples", "qOtherOut" )
> 
> # mean-scaling input and output quantities
> for( n in c( xNames, yNames ) ) {
+   appleProdFr86[[ n ]] <- appleProdFr86[[ n ]] / mean( appleProdFr86[[ n ]] ) 
+ }
> 
> 
> ### Cobb-Douglas ray-based input distance function
> estCD <- distRayEst( xNames = xNames, yNames = yNames,
+   data = appleProdFr86, form = "cd" )
Warning message:
The residuals of the OLS are right-skewed. This may indicate the absence of inefficiency or
  model misspecification or sample 'bad luck' 
> cbind( round( coef( estCD )[ !grepl( "(Intercept)", names( coef( estCD ) ) ) ], 2 ) )
         [,1]
alpha_1  0.09
alpha_2  0.62
beta_1   0.17
beta_2  -0.39
> ## IGNORE_RDIFF_BEGIN
> cbind( round( coef( estCD )[ grepl( "(Intercept)", names( coef( estCD ) ) ) ], 2 ) )
                [,1]
(Intercept)     0.07
Zu_(Intercept) -9.83
Zv_(Intercept) -2.20
> cbind( names( estCD ) )
      [,1]          
 [1,] "call"        
 [2,] "formula"     
 [3,] "S"           
 [4,] "typeSfa"     
 [5,] "Nobs"        
 [6,] "nXvar"       
 [7,] "scaling"     
 [8,] "logDepVar"   
 [9,] "nuZUvar"     
[10,] "nvZVvar"     
[11,] "nParm"       
[12,] "udist"       
[13,] "startVal"    
[14,] "dataTable"   
[15,] "olsParam"    
[16,] "olsStder"    
[17,] "olsSigmasq"  
[18,] "olsLoglik"   
[19,] "olsSkew"     
[20,] "olsM3Okay"   
[21,] "CoelliM3Test"
[22,] "AgostinoTest"
[23,] "isWeights"   
[24,] "optType"     
[25,] "nIter"       
[26,] "optStatus"   
[27,] "startLoglik" 
[28,] "mlLoglik"    
[29,] "mlParam"     
[30,] "gradient"    
[31,] "gradL_OBS"   
[32,] "gradientNorm"
[33,] "invHessian"  
[34,] "hessianType" 
[35,] "mlDate"      
[36,] "coefList"    
[37,] "ela"         
[38,] "mono"        
> ## IGNORE_RDIFF_END
> lapply( estCD$coefList, function(x) round( x, 2 ) )
$alpha0
(Intercept) 
       0.07 

$alphaVec
[1] 0.09 0.62 0.30

$betaVec
[1]  0.17 -0.39

> apply( estCD$mono, 2, table )
$x_1

TRUE 
 140 

$x_2

TRUE 
 140 

$x_3

TRUE 
 140 

$y_1

TRUE 
 140 

$y_2

FALSE  TRUE 
   25   115 

$all_x

TRUE 
 140 

$all_y

FALSE  TRUE 
   25   115 

$all

FALSE  TRUE 
   25   115 

> lapply( estCD$ela, function(x) round( summary(x), 2 ) )
$ela_x_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.09    0.09    0.09    0.09    0.09    0.09 

$ela_x_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.62    0.62    0.62    0.62    0.62    0.62 

$ela_x_3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.3     0.3     0.3     0.3     0.3     0.3 

$ela_y_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.41   -0.33   -0.14   -0.18   -0.04    0.00 

$ela_y_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.39   -0.35   -0.25   -0.21   -0.06    0.02 

$ela_scale
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.56    2.56    2.56    2.56    2.56    2.56 

> 
> # calculate the dependent variable (logarithm of predicted distance)
> appleProdFr86$distCD <- distRayCalc( xNames = xNames, yNames = yNames, 
+   data = appleProdFr86, coef = coef( estCD ), form = "cd" )
> round( summary( appleProdFr86$distCD ), 2 )
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.76   -0.23    0.01    0.01    0.24    0.76 
> 
> # calculate elasticities
> elaCD <- distRayEla( xNames = xNames, yNames = yNames,
+   coef = coef( estCD ), data = appleProdFr86, form = "cd" )
> all.equal( elaCD, estCD$ela )
[1] TRUE
> lapply( elaCD, function(x) round( summary(x), 2 ) )
$ela_x_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.09    0.09    0.09    0.09    0.09    0.09 

$ela_x_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.62    0.62    0.62    0.62    0.62    0.62 

$ela_x_3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.3     0.3     0.3     0.3     0.3     0.3 

$ela_y_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.41   -0.33   -0.14   -0.18   -0.04    0.00 

$ela_y_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.39   -0.35   -0.25   -0.21   -0.06    0.02 

$ela_scale
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.56    2.56    2.56    2.56    2.56    2.56 

> 
> # calculate derivatives
> derivCD <- distRayDeriv( xNames = xNames, yNames = yNames,
+   coef = coef( estCD ), data = appleProdFr86, form = "cd" )
> lapply( derivCD, function(x) round( summary(x), 2 ) )
$d_x_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.02    0.07    0.11    0.16    0.17    1.20 

$d_x_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.24    0.56    0.73    0.78    0.91    2.09 

$d_x_3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.12    0.26    0.35    0.44    0.50    1.78 

$d_y_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -6.07   -0.44   -0.27   -0.46   -0.15   -0.02 

$d_y_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -6.82   -0.40   -0.19   -0.44   -0.05    0.15 

> 
> # vector of unrestricted coefficients and their covariance matrix
> nCoefCD <- length( coef( estCD ) ) - 2
> uCoefCD <- coef( estCD )[ 1:nCoefCD ]
> uCovInvCD <- solve( vcov( estCD )[ 1:nCoefCD, 1:nCoefCD ] )
> 
> # obtain the matrix and vector to impose monotonicity
> restrCD <- distRayMonoRestr( xNames = xNames, yNames = yNames, 
+   data = appleProdFr86, form = "cd" )
> 
> # obtain the restricted coefficients
> minDistCD <- solve.QP( Dmat = uCovInvCD, dvec = rep( 0, nCoefCD ),
+   Amat = t( restrCD$RMat ), bvec = - restrCD$RMat %*% uCoefCD + restrCD$rVec )
> rCoefCD <- minDistCD$solution + uCoefCD
> round( rCoefCD, 2 )
(Intercept)     alpha_1     alpha_2      beta_1      beta_2 
       0.23        0.09        0.62        0.00       -0.42 
> 
> # calculate elasticities based on restricted coefficients
> rElaCD <- distRayEla( xNames = xNames, yNames = yNames,
+   coef = rCoefCD, data = appleProdFr86, form = "cd" )
> lapply( rElaCD, function(x) round( summary(x), 2 ) )
$ela_x_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.09    0.09    0.09    0.09    0.09    0.09 

$ela_x_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.62    0.62    0.62    0.62    0.62    0.62 

$ela_x_3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.29    0.29    0.29    0.29    0.29    0.29 

$ela_y_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.42   -0.27   -0.08   -0.15   -0.01    0.00 

$ela_y_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.42   -0.41   -0.34   -0.27   -0.15    0.00 

$ela_scale
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.38    2.38    2.38    2.38    2.38    2.38 

> 
> 
> ### Translog ray-based input distance function
> estTL <- distRayEst( xNames = xNames, yNames = yNames,
+   data = appleProdFr86 )
> cbind( round( coef( estTL ), 2 ) )
                [,1]
(Intercept)     0.83
alpha_1         0.11
alpha_2         0.44
beta_1         -0.94
beta_2         -0.77
alpha_1_1       0.16
alpha_1_2       0.03
alpha_2_2       0.04
beta_1_1        1.10
beta_1_2        0.37
beta_2_2       -0.06
psi_1_1        -0.02
psi_1_2        -0.02
psi_2_1         0.21
psi_2_2         0.06
Zu_(Intercept) -2.01
Zv_(Intercept) -3.47
> ## IGNORE_RDIFF_BEGIN
> cbind( names( estTL ) )
      [,1]          
 [1,] "call"        
 [2,] "formula"     
 [3,] "S"           
 [4,] "typeSfa"     
 [5,] "Nobs"        
 [6,] "nXvar"       
 [7,] "scaling"     
 [8,] "logDepVar"   
 [9,] "nuZUvar"     
[10,] "nvZVvar"     
[11,] "nParm"       
[12,] "udist"       
[13,] "startVal"    
[14,] "dataTable"   
[15,] "olsParam"    
[16,] "olsStder"    
[17,] "olsSigmasq"  
[18,] "olsLoglik"   
[19,] "olsSkew"     
[20,] "olsM3Okay"   
[21,] "CoelliM3Test"
[22,] "AgostinoTest"
[23,] "isWeights"   
[24,] "optType"     
[25,] "nIter"       
[26,] "optStatus"   
[27,] "startLoglik" 
[28,] "mlLoglik"    
[29,] "mlParam"     
[30,] "gradient"    
[31,] "gradL_OBS"   
[32,] "gradientNorm"
[33,] "invHessian"  
[34,] "hessianType" 
[35,] "mlDate"      
[36,] "coefList"    
[37,] "ela"         
[38,] "mono"        
> ## IGNORE_RDIFF_END
> lapply( estTL$coefList, function(x) round( x, 2 ) )
$alpha0
(Intercept) 
       0.83 

$alphaVec
[1] 0.11 0.44 0.45

$betaVec
[1] -0.94 -0.77

$alphaMat
      [,1]  [,2]  [,3]
[1,]  0.16  0.03 -0.19
[2,]  0.03  0.04 -0.07
[3,] -0.19 -0.07  0.26

$betaMat
     [,1]  [,2]
[1,] 1.10  0.37
[2,] 0.37 -0.06

$psiMat
      [,1]  [,2]
[1,] -0.02 -0.02
[2,]  0.21  0.06
[3,] -0.19 -0.04

> apply( estTL$mono, 2, table )
$x_1

FALSE  TRUE 
   36   104 

$x_2

TRUE 
 140 

$x_3

FALSE  TRUE 
    8   132 

$y_1

FALSE  TRUE 
    6   134 

$y_2

FALSE  TRUE 
    1   139 

$all_x

FALSE  TRUE 
   44    96 

$all_y

FALSE  TRUE 
    7   133 

$all

FALSE  TRUE 
   47    93 

> lapply( estTL$ela, function(x) round( summary(x), 2 ) )
$ela_x_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.19    0.00    0.07    0.07    0.14    0.31 

$ela_x_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.40    0.58    0.66    0.66    0.74    0.88 

$ela_x_3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.12    0.15    0.28    0.27    0.37    0.69 

$ela_y_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.91   -0.36   -0.16   -0.24   -0.06    0.15 

$ela_y_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.60   -0.24   -0.16   -0.17   -0.10    0.01 

$ela_scale
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.09    1.87    2.75    4.08    3.95  112.32 

> 
> # calculate the dependent variable (logarithm of predicted distance)
> appleProdFr86$logDistTL <- distRayCalc( xNames = xNames, yNames = yNames, 
+   data = appleProdFr86, coef = coef( estTL ) )
> round( summary( appleProdFr86$logDistTL ), 2 )
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.40    0.09    0.26    0.29    0.47    1.09 
> 
> # calculate elasticities
> elaTL <- distRayEla( xNames = xNames, yNames = yNames,
+   coef = coef( estTL ), data = appleProdFr86 )
> all.equal( elaTL, estTL$ela )
[1] TRUE
> lapply( elaTL, function(x) round( summary(x), 2 ) )
$ela_x_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.19    0.00    0.07    0.07    0.14    0.31 

$ela_x_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.40    0.58    0.66    0.66    0.74    0.88 

$ela_x_3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.12    0.15    0.28    0.27    0.37    0.69 

$ela_y_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.91   -0.36   -0.16   -0.24   -0.06    0.15 

$ela_y_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.60   -0.24   -0.16   -0.17   -0.10    0.01 

$ela_scale
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.09    1.87    2.75    4.08    3.95  112.32 

> 
> # calculate derivatives
> derivTL <- distRayDeriv( xNames = xNames, yNames = yNames,
+   coef = coef( estTL ), data = appleProdFr86 )
> lapply( derivTL, function(x) round( summary(x), 
+   3 - round( log( max( abs( x ) ), 10 ) ) ) )
$d_x_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -4.52    0.00    0.10    0.01    0.21    0.89 

$d_x_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.166   0.887   1.134   1.110   1.399   2.275 

$d_x_3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -0.508   0.229   0.389   0.457   0.617   2.630 

$d_y_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -1.82   -0.90   -0.48   -0.51   -0.30    7.23 

$d_y_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -9.33   -0.60   -0.27   -0.59   -0.13    0.01 

> 
> # vector of unrestricted coefficients and their covariance matrix
> nCoefTL <- length( coef( estTL ) ) - 2
> uCoefTL <- coef( estTL )[ 1:nCoefTL ]
> uCovInvTL <- solve( vcov( estTL )[ 1:nCoefTL, 1:nCoefTL ] )
> 
> # obtain the matrix and vector to impose monotonicity
> restrTL <- distRayMonoRestr( xNames = xNames, yNames = yNames, 
+   data = appleProdFr86 )
> 
> # obtain the restricted coefficientslibrary( "quadprog" )
> minDistTL <- solve.QP( Dmat = uCovInvTL, dvec = rep( 0, nCoefTL ),
+   Amat = t( restrTL$RMat ), bvec = - restrTL$RMat %*% uCoefTL + restrTL$rVec )
> rCoefTL <- minDistTL$solution + uCoefTL
> round( rCoefTL, 2 )
(Intercept)     alpha_1     alpha_2      beta_1      beta_2   alpha_1_1 
       0.73        0.07        0.44       -0.77       -0.65        0.04 
  alpha_1_2   alpha_2_2    beta_1_1    beta_1_2    beta_2_2     psi_1_1 
       0.03        0.01        1.01        0.24       -0.07        0.01 
    psi_1_2     psi_2_1     psi_2_2 
       0.01        0.20        0.08 
> 
> # calculate elasticities based on restricted coefficients
> rElaTL <- distRayEla( xNames = xNames, yNames = yNames,
+   coef = rCoefTL, data = appleProdFr86 )
> lapply( rElaTL, function(x) round( summary(x), 2 ) )
$ela_x_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    0.06    0.09    0.08    0.11    0.16 

$ela_x_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.39    0.57    0.63    0.64    0.71    0.85 

$ela_x_3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    0.20    0.29    0.28    0.36    0.61 

$ela_y_1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.79   -0.35   -0.19   -0.24   -0.10    0.00 

$ela_y_2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.45   -0.25   -0.18   -0.18   -0.11    0.00 

$ela_scale
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.26    1.86    2.52    2.94    3.23   17.05 

> 
> 
> proc.time()
   user  system elapsed 
   2.00    0.15    2.15