Data importation and variable definitions

library("readxl")
my_data <- read_excel("mydatax.xlsx")
m = my_data$Markup
mm = seq(from = min(m), to = max(m), length.out = 10^4)
mmin = min(m)
mn = max(m)
nn = length(m)

Table 2: Estimated Markup Densities Given Assumptions About Productivity and Demand

Pareto productivity distribution + CREMR demand

# Objective function (i.e.  minus log-likelihood) to be minimized:
log_likelihood_b_CREMR = function(theta, m){
  k = theta[1]
  sigma = theta[2]
  -sum(log((k*(mmin/(mmin+sigma-sigma*mmin))^(k/sigma)/(m*(m+sigma-m*sigma)))*(m/(m+sigma-m*sigma))^(-k/sigma)))
}
# Estimation (CREMR+Pareto)
theta_startCREMR = c(2, 0.5*(1 + mn/(mn-1)))
(estimCREMR = optim(par = theta_startCREMR, log_likelihood_b_CREMR, m = m))
## $par
## [1] 1.232891 1.111174
## 
## $value
## [1] 3021.716
## 
## $counts
## function gradient 
##       87       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
(kCREMR=estimCREMR$par[1])
## [1] 1.232891
(sigmaCREMR=estimCREMR$par[2])
## [1] 1.111174
(omegaCREMR=((sigmaCREMR-1)*mmin/(mmin+sigmaCREMR-sigmaCREMR*mmin))^(1/sigmaCREMR))
## [1] 0.1386922
# AIC (with three estimated parameters)
pCREMR = 3
(AIC_CREMR = 2*pCREMR - 2*(-estimCREMR$value))
## [1] 6049.433

Pareto productivity distribution + Linear demand

# MLE estimators (closed form)
(omegaLIN=2*mmin-1)
## [1] 1.002774
(kkLIN = nn/(sum(log(2*m-1))-nn*log(omegaLIN)))
## [1] 1.001127
# AIC (with two estimated parameters):
pLIN = 2
(AIC_LIN = 2*pLIN - 2*(sum(log(2*kkLIN*(omegaLIN^kkLIN)*(2*m-1)^(-kkLIN-1)))))
## [1] 6428.425

Pareto productivity distribution + LES demand

# MLE estimators (closed form)
(omegaLES=mmin^2)
## [1] 1.002776
(kkLES = 1/((2/nn)*sum(log(m))-log(omegaLES)))
## [1] 0.7466883
# AIC (with two estimated parameters):
pLES = 2
(AIC_LES = 2*pLES - 2*(sum(log((2*kkLES*(omegaLES^(kkLES))*(m)^(-2*kkLES-1))))))
## [1] 6244.631

Pareto productivity distribution + Translog demand

# MLE estimators (closed form)
(omegaTLOG=mmin*exp(mmin))
## [1] 2.72583
(kkTLOG = 1/((1/nn)*sum(log(m)+m)-log(omegaTLOG)))
## [1] 0.4868533
# AIC (with two estimated parameters):
pTLOG = 2
(AIC_TLOG = 2*pTLOG - 2*(sum(log(kkTLOG*(omegaTLOG^(kkTLOG))*(1+m)*m^(-kkTLOG-1)*exp(-kkTLOG*m)))))
## [1] 6258.079

Figure 4(a): Comparing the four demand alternatives (Pareto productivity)

hist(m, probability = TRUE, col = "lightgrey", breaks = 20)
mm = seq(from = min(m), to = max(m), length.out = 10^4)
xx = 2*kkLIN*(omegaLIN^kkLIN)*(2*mm-1)^(-kkLIN-1)
lines(mm, xx, col = 3, lwd = 2)
vv = 2*kkLES*omegaLES^(kkLES)*(mm)^(-2*kkLES-1)
lines(mm, vv, col = 4, lwd = 2)
ww = kkTLOG*(omegaTLOG^(kkTLOG))*(1+mm)*mm^(-kkTLOG-1)*exp(-kkTLOG*mm)
lines(mm, ww, col = 5, lwd = 2)
yy = (kCREMR/mm^2)*(mm/(mm+sigmaCREMR-mm*sigmaCREMR))^((sigmaCREMR-kCREMR)/sigmaCREMR)
lines(mm, yy, col = 2, lwd = 2)
legend("topright", c("CREMR","Linear","LES","Translog"), col = 2:5, lwd = 2)

Lognormal productivity distribution + CREMR demand

# Auxiliary function to ensure estimation over the appropriate range
a=mn/(mn-1)
u=function(z){
  ((a-1)/pi)*atan(z)+(a+1)/2
}
# Objective function (i.e.  minus log-likelihood) to be minimized (untruncated lognormal)
log_likelihood_b_CREMR_ln = function(theta, m){
  mu = theta[1]
  z = theta[2]
  ss=theta[3]
  -sum(log((exp(-((log((m*(u(z)-1))/(m+u(z)-m*u(z)))-mu)^2/(2*(exp(ss)*u(z))^2))))/(sqrt(2*pi)*m*exp(ss)*(m+u(z)-m*u(z)))))
}
# Estimation (untruncated lognormal)
# Starting values
sigma_start=0.5*(1 + mn/(mn-1))
z_start=tan((pi/(a-1))*(sigma_start-(a+1)/2))
theta_start = c(2, z_start, 1)
(estimCREMR_ln = optim(par = theta_start, log_likelihood_b_CREMR_ln,method = "BFGS",m = m))
## $par
## [1]  -6.1145127 -31.4177698  -0.5466739
## 
## $value
## [1] 3788.063
## 
## $counts
## function gradient 
##      365       97 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
# Estimated parameters (untruncated lognormal)
(muCREMR_ln=estimCREMR_ln$par[1])
## [1] -6.114513
(zCREMR_ln=estimCREMR_ln$par[2])
## [1] -31.41777
(ssCREMR_ln=estimCREMR_ln$par[3])
## [1] -0.5466739
(sigmaCREMR_ln=u(zCREMR_ln))
## [1] 1.001128
(sCREMR_ln=exp(ssCREMR_ln))
## [1] 0.578872
# Truncated lognormal case
# Truncation point
(mt=min(m))
## [1] 1.001387
# Objective function (i.e.  minus log-likelihood) to be minimized (truncated lognormal)
log_likelihood_b_CREMR_lnt = function(theta, m){
  mu = theta[1]
  z = theta[2]
  s=theta[3]
  -sum(log((1/(1-pnorm((log((mt*(u(z)-1))/(mt+u(z)-mt*u(z)))-mu)/(u(z)*exp(s)),0,1)))*(exp(-((log((m*(u(z)-1))/(m+u(z)-m*u(z)))-mu)^2/(2*(exp(s)*u(z))^2))))/(sqrt(2*pi)*m*exp(s)*(m+u(z)-m*u(z)))))
}
# Estimation (truncated lognormal)
# Starting values from untruncated estimation
theta_start_t = c(muCREMR_ln, zCREMR_ln, ssCREMR_ln)
(estimCREMR_lnt = optim(par = theta_start, log_likelihood_b_CREMR_lnt,m = m))
## $par
## [1] -49.982070  29.043672   1.800142
## 
## $value
## [1] 3026.494
## 
## $counts
## function gradient 
##      347       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
# Estimated parameters (truncated lognormal)
(muCREMR_lnt=estimCREMR_lnt$par[1])
## [1] -49.98207
(zCREMR_lnt=estimCREMR_lnt$par[2])
## [1] 29.04367
(ssCREMR_lnt=estimCREMR_lnt$par[3])
## [1] 1.800142
(sigmaCREMR_lnt=u(zCREMR_lnt))
## [1] 1.110173
(sCREMR_lnt=exp(ssCREMR_lnt))
## [1] 6.050505
(omegaCREMR_lnt=((mt*(sigmaCREMR_lnt-1))/(mt+sigmaCREMR_lnt-mt*sigmaCREMR_lnt))^(1/sigmaCREMR_lnt))
## [1] 0.1373217
(truncCREMR_lnt=pnorm((log((mt*(sigmaCREMR_lnt-1))/(mt+sigmaCREMR_lnt-mt*sigmaCREMR_lnt))-muCREMR_lnt)/(sigmaCREMR_lnt*sCREMR_lnt),0,1))
## [1] 1
(1/(1-pnorm((log((mt*(sigmaCREMR_lnt-1))/(mt+sigmaCREMR_lnt-mt*sigmaCREMR_lnt))-muCREMR_lnt)/(sigmaCREMR_lnt*sCREMR_lnt),0,1)))
## [1] 1.759906e+12
(1/(1-truncCREMR_lnt))
## [1] 1.759906e+12
(formatC(1/(1-truncCREMR_lnt), digits = 4, format = "f"))
## [1] "1759906067749.3147"
# AIC (with four estimated parameters):
pCREMR_lnt = 4
(AIC_CREMR_lnt = 2*pCREMR_lnt - 2*(-estimCREMR_lnt$value))
## [1] 6060.988

Lognormal productivity distribution + Linear demand

# MLE estimation (untruncated lognormal)
(muLIN=(1/nn)*sum(log(2*m-1)))
## [1] 1.001645
(sLIN=sqrt((1/nn)*sum((log(2*m-1)-muLIN)^2)))
## [1] 0.7497347
# MLE estimation (truncated lognormal)
log_likelihood_h_LIN_lnt = function(theta, m){
  muLINt = theta[1]
  sLINt = theta[2]
  -sum(log((1/(1-pnorm((log(2*min(m)-1)-muLINt)/sLINt, 0, 1)))*sqrt(2/pi)*(1/(sLINt*(2*m-1)))*exp(-((log(2*m-1)-muLINt)^2)/(2*sLINt^2))))
}
# Starting values from untruncated estimation
theta_startLIN_lnt = c(muLIN, sLIN)
(estimLIN_lnt = optim(par = theta_startLIN_lnt, log_likelihood_h_LIN_lnt, m = m))
## $par
## [1] 0.05394534 1.22821098
## 
## $value
## [1] 3087.331
## 
## $counts
## function gradient 
##       51       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
# AIC (with three estimated parameters):
pLIN_lnt = 3
(AIC_LIN_lnt = 2*pLIN_lnt - 2*(-estimLIN_lnt$value))
## [1] 6180.661

Lognormal productivity distribution + LES demand

# MLE estimation (untruncated lognormal)
(muLES=(2/nn)*sum(log(m)))
## [1] 1.342019
(sLES=sqrt((1/nn)*sum((2*log(m)-muLES)^2)))
## [1] 1.155128
# MLE estimation (truncated lognormal)
log_likelihood_h_LES_lnt = function(theta, m){
  muLESt = theta[1]
  sLESt = theta[2]
  -sum(log((1/(1-pnorm((2*log(min(m))-muLESt)/sLESt, 0, 1)))*sqrt(2/pi)*(1/(sLESt*m))*exp(-((2*log(m)-muLESt)^2)/(2*sLESt^2))))
}
# Starting values from untruncated estimation
theta_startLES_lnt = c(muLES, sLES)
(estimLES_lnt = optim(par = theta_startLES_lnt, log_likelihood_h_LES_lnt, m = m))
## $par
## [1] -3.234377  2.731872
## 
## $value
## [1] 3089.392
## 
## $counts
## function gradient 
##       75       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
# AIC (with three estimated parameters):
pLES_lnt = 3
(AIC_LES_lnt = 2*pLES_lnt - 2*(-estimLES_lnt$value))
## [1] 6184.785

Lognormal productivity distribution + Translog demand

# MLE estimation (untruncated lognormal) - closed form
(muTLOG=(1/nn)*sum(log(m*exp(m))))
## [1] 3.05678
(sTLOG=sqrt((1/nn)*sum((log(m*exp(m))-muTLOG)^2)))
## [1] 2.39093
# MLE estimation (truncated lognormal)
log_likelihood_h_TLOG_lnt = function(theta, m){
  mu = theta[1]
  s = theta[2]
  mt=min(m)
  -sum(log((1/(1-pnorm((log(mt*exp(mt))-mu)/exp(s), 0, 1)))*(1/(sqrt(2*pi)*exp(s)))*((m+1)/m)*exp(-((log(m*exp(m))-mu)^2)/(2*exp(s)^2))))
}
# Starting values from untruncated estimation
theta_startTLOG_lnt = c(muTLOG, log(sTLOG))
(estimTLOG_lnt = optim(par = theta_startTLOG_lnt, log_likelihood_h_TLOG_lnt,method = "BFGS",m = m))
## $par
## [1] -77.641030   2.569909
## 
## $value
## [1] 3138.547
## 
## $counts
## function gradient 
##      167       46 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
# Estimated parameters (truncated lognormal)
(muTLOGt=estimTLOG_lnt$par[1])
## [1] -77.64103
(sTLOGt=exp(estimTLOG_lnt$par[2]))
## [1] 13.06463
# AIC (with three estimated parameters)
pTLOG_lnt = 3
(AIC_TLOG_lnt = 2*pTLOG_lnt - 2*(-estimTLOG_lnt$value))
## [1] 6283.095

Figure 4(b): Comparing the four demand alternatives (truncated lognormal productivity)

hist(m, probability = TRUE, col = "lightgrey", breaks = 20)
mm = seq(from = min(m), to = max(m), length.out = 10^4)
ww = (1/(1-pnorm((log(mt*exp(mt))-muTLOGt)/sTLOGt, 0, 1)))*(1/(sqrt(2*pi)*sTLOGt))*((mm+1)/mm)*exp(-((log(mm*exp(mm))-muTLOGt)^2)/(2*sTLOGt^2))
lines(mm, ww, col = 5, lwd = 2)
xx = (1/(1-pnorm((log(2*min(m)-1)-estimLIN_lnt$par[1])/estimLIN_lnt$par[2], 0, 1)))*sqrt(2/pi)*(1/(estimLIN_lnt$par[2]*(2*mm-1)))*exp(-((log(2*mm-1)-estimLIN_lnt$par[1])^2)/(2*estimLIN_lnt$par[2]^2))
lines(mm, xx, col = 3, lwd = 2)
vv = (1/(1-pnorm((2*log(min(m))-estimLES_lnt$par[1])/estimLES_lnt$par[2], 0, 1)))*sqrt(2/pi)*(1/(estimLES_lnt$par[2]*mm))*exp(-((2*log(mm)-estimLES_lnt$par[1])^2)/(2*estimLES_lnt$par[2]^2))
lines(mm, vv, col = 4, lwd = 2)
yy = ((1/(1-pnorm((log((mt*(sigmaCREMR_lnt-1))/(mt+sigmaCREMR_lnt-mt*sigmaCREMR_lnt))-muCREMR_lnt)/(sigmaCREMR_lnt*sCREMR_lnt),0,1))))*(exp(-((log((mm*sigmaCREMR_lnt-mm)/(sigmaCREMR_lnt+mm-sigmaCREMR_lnt*mm))-muCREMR_lnt)^2/(2*(sCREMR_lnt*sigmaCREMR_lnt)^2))))/(sqrt(2*pi)*mm*sCREMR_lnt*(mm+sigmaCREMR_lnt-mm*sigmaCREMR_lnt))
lines(mm, yy, col = 2, lwd = 2)
legend("topright", c("CREMR","Linear","LES","Translog"), col = 2:5, lwd = 2)

AIC ranking for all 8 cases

library(kableExtra)
dt = data.frame(models = c("Pareto - CREMR", "LN - CREMR","LN - Linear", "LN - LES", "Pareto-LES","Pareto-Translog", "LN - Translog","Pareto-Linear"), AIC = c(AIC_CREMR,AIC_CREMR_lnt,AIC_LIN_lnt, AIC_LES_lnt,AIC_LES,AIC_TLOG,AIC_TLOG_lnt,AIC_LIN))
kable(dt)
models AIC
Pareto - CREMR 6049.433
LN - CREMR 6060.988
LN - Linear 6180.661
LN - LES 6184.785
Pareto-LES 6244.631
Pareto-Translog 6258.079
LN - Translog 6283.095
Pareto-Linear 6428.425

Misallocation calculations

Underlying productivity distribution: Pareto

CDF and PDF definitions

CDFxMarketP = function(sigma,gama,omega,k, x){
  1-gama^(k/sigma)*omega^k*(x-gama)^(-k/sigma)
}
PDFxMarketP=function(sigma,gama,omega,k, x){
  gama^(k/sigma)*(k/sigma)*omega^k*(x-gama)^(-(k+sigma)/sigma)
  }
CDFxPlanP = function(sigma,gama,omega,k, x){
  1-gama^(k/sigma)*((omega^(sigma-1)/(1+omega^sigma))*x*(x-gama)^((1-sigma)/sigma))^(-k)
}
PDFxPlanP = function(sigma,gama,omega,k, x){
 gama^(k/sigma)*(k/sigma)*((x-gama*sigma)/(x*(x-gama)))*((omega^(sigma-1)/(1+omega^sigma))*x*(x-gama)^((1-sigma)/sigma))^(-k)
}
PDFdifferenceP=function(sigma,gama,omega,k, x){
  PDFxPlanP(sigma,gama,omega,k, x)-PDFxMarketP(sigma,gama,omega,k, x)
}

Calculating the intersections between PDFs and shares:

gama1=1
f1P=function(x){
  PDFdifferenceP(sigmaCREMR,gama1,omegaCREMR,kCREMR, x)
}
uniroot(f1P,c(1.4,1.5))
## $root
## [1] 1.464798
## 
## $f.root
## [1] 3.223167e-05
## 
## $iter
## [1] 4
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc1P<-uniroot(f1P,c(1.4,1.5),tol = 1e-6)
CDFxMarketP(sigmaCREMR,gama1,omegaCREMR,kCREMR, xc1P$root)
## [1] 0.7951459
CDFxPlanP(sigmaCREMR,gama1,omegaCREMR,kCREMR, xc1P$root)
## [1] 0.1514132
CDFxMarketP(sigmaCREMR,gama1,omegaCREMR,kCREMR, xc1P$root)/CDFxPlanP(sigmaCREMR,gama1,omegaCREMR,kCREMR, xc1P$root)
## [1] 5.251495
gama3=1.5
f3P=function(x){
  PDFdifferenceP(sigmaCREMR,gama3,omegaCREMR,kCREMR, x)
}
uniroot(f3P,c(2.1,2.2))
## $root
## [1] 2.19717
## 
## $f.root
## [1] -1.460128e-06
## 
## $iter
## [1] 3
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc3P<-uniroot(f3P,c(2.1,2.2),tol = 1e-6)
CDFxMarketP(sigmaCREMR,gama3,omegaCREMR,kCREMR, xc3P$root)
## [1] 0.7951459
CDFxPlanP(sigmaCREMR,gama3,omegaCREMR,kCREMR, xc3P$root)
## [1] 0.1514132
CDFxMarketP(sigmaCREMR,gama3,omegaCREMR,kCREMR, xc3P$root)/CDFxPlanP(sigmaCREMR,gama3,omegaCREMR,kCREMR, xc3P$root)
## [1] 5.251496
gama5=2
f5P=function(x){
  PDFdifferenceP(sigmaCREMR,gama5,omegaCREMR,kCREMR, x)
}
uniroot(f5P,c(2.9,2.95))
## $root
## [1] 2.929558
## 
## $f.root
## [1] -1.993847e-06
## 
## $iter
## [1] 3
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc5P<-uniroot(f5P,c(2.9,2.95),tol = 1e-6)
CDFxMarketP(sigmaCREMR,gama5,omegaCREMR,kCREMR, xc5P$root)
## [1] 0.7951459
CDFxPlanP(sigmaCREMR,gama5,omegaCREMR,kCREMR, xc5P$root)
## [1] 0.1514132
CDFxMarketP(sigmaCREMR,gama5,omegaCREMR,kCREMR, xc5P$root)/CDFxPlanP(sigmaCREMR,gama5,omegaCREMR,kCREMR, xc5P$root)
## [1] 5.251496
gama6=3
f6P=function(x){
  PDFdifferenceP(sigmaCREMR,gama6,omegaCREMR,kCREMR, x)
}
uniroot(f6P,c(4.38,4.4))
## $root
## [1] 4.394343
## 
## $f.root
## [1] -4.759014e-09
## 
## $iter
## [1] 2
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 0.0001121995
xc6P<-uniroot(f6P,c(4.38,4.4),tol = 1e-6)
CDFxMarketP(sigmaCREMR,gama6,omegaCREMR,kCREMR, xc6P$root)
## [1] 0.7951459
CDFxPlanP(sigmaCREMR,gama6,omegaCREMR,kCREMR, xc6P$root)
## [1] 0.1514132
CDFxMarketP(sigmaCREMR,gama6,omegaCREMR,kCREMR, xc6P$root)/CDFxPlanP(sigmaCREMR,gama6,omegaCREMR,kCREMR, xc6P$root)
## [1] 5.251496
gama7=4
f7P=function(x){
  PDFdifferenceP(sigmaCREMR,gama7,omegaCREMR,kCREMR, x)
}
uniroot(f7P,c(5.84,5.86))
## $root
## [1] 5.859141
## 
## $f.root
## [1] 2.097688e-06
## 
## $iter
## [1] 2
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc7P<-uniroot(f7P,c(5.84,5.86),tol = 1e-6)
CDFxMarketP(sigmaCREMR,gama7,omegaCREMR,kCREMR, xc7P$root)
## [1] 0.7951459
CDFxPlanP(sigmaCREMR,gama7,omegaCREMR,kCREMR, xc7P$root)
## [1] 0.1514132
CDFxMarketP(sigmaCREMR,gama7,omegaCREMR,kCREMR, xc7P$root)/CDFxPlanP(sigmaCREMR,gama7,omegaCREMR,kCREMR, xc7P$root)
## [1] 5.251496
gama8=5
f8P=function(x){
  PDFdifferenceP(sigmaCREMR,gama8,omegaCREMR,kCREMR, x)
}
uniroot(f8P,c(7.31,7.33))
## $root
## [1] 7.323904
## 
## $f.root
## [1] -4.053606e-10
## 
## $iter
## [1] 2
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 7.041729e-05
xc8P<-uniroot(f8P,c(7.31,7.33),tol = 1e-6)
CDFxMarketP(sigmaCREMR,gama8,omegaCREMR,kCREMR, xc8P$root)
## [1] 0.7951459
CDFxPlanP(sigmaCREMR,gama8,omegaCREMR,kCREMR, xc8P$root)
## [1] 0.1514132
CDFxMarketP(sigmaCREMR,gama8,omegaCREMR,kCREMR, xc8P$root)/CDFxPlanP(sigmaCREMR,gama8,omegaCREMR,kCREMR, xc8P$root)
## [1] 5.251496

Underlying productivity distribution: truncated lognormal

CDF and PDF definitions

CDFxMarketLN = function(sigma,gama,mu,s, x){
  pnorm((log(x-gama)-mu)/(sigma*s),0,1)
}
CDFxMarketTLN = function(sigma,gama,mu,s,xbar,x){
  (pnorm((log(x-gama)-mu)/(sigma*s),0,1)-pnorm((log(xbar-gama)-mu)/(sigma*s),0,1))/(1-pnorm((log(xbar-gama)-mu)/(sigma*s),0,1))
}
PDFxMarketTLN=function(sigma,gama,mu,s,xbar,x){
  (1/(1-pnorm((log(xbar-gama)-mu)/(sigma*s),0,1)))*(exp(-(mu-log(x-gama))^2/(2*s^2*sigma^2))/(sqrt(2*pi)*s*(x-gama)*sigma))
  }
CDFxPlanLN = function(sigma,gama,mu,s,omega,x){
  pnorm((log((omega^sigma/(1+omega^sigma))*x*(x-gama)^((1-sigma)/sigma))-mu)/(s),0,1)
}
CDFxPlanTLN = function(sigma,gama,mu,s,omega,xbar,x){
  (pnorm((log((omega^sigma/(1+omega^sigma))*x*(x-gama)^((1-sigma)/sigma))-mu)/(s),0,1)-pnorm((log((omega^sigma/(1+omega^sigma))*xbar*(xbar-gama)^((1-sigma)/sigma))-mu)/(s),0,1))/(1-pnorm((log((omega^sigma/(1+omega^sigma))*xbar*(xbar-gama)^((1-sigma)/sigma))-mu)/(s),0,1))
}
PDFxPlanTLN = function(sigma,gama,mu,s,omega,xbar, x){
  (1/(1-pnorm(((log((omega^sigma/(1+omega^sigma))*xbar*(xbar-gama)^((1-sigma)/sigma))-mu)/s),0,1)))*(exp(-(mu-log((omega^sigma/(1+omega^sigma))*x*(x-gama)^(1/sigma-1)))^2/(2*s^2))*(x-gama*sigma)/(sqrt(2*pi)*s*x*(x-gama)*sigma))
}
PDFdifferenceLN=function(sigma,gama,muM,muP,s,omega,xbar,x){
  PDFxPlanTLN(sigma,gama,muP,s,omega,xbar, x)-PDFxMarketTLN(sigma,gama,muM,s,xbar, x)
}

Calculating the intersections between PDFs and shares:

gama1=1
xbar1=gama1*(omegaCREMR_lnt^sigmaCREMR_lnt+1)
f1LN=function(x){
  PDFdifferenceLN(sigmaCREMR_lnt,gama1,muCREMR_lnt+log(gama1),(muCREMR_lnt+log(gama1))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar1,x)
}
uniroot(f1LN,c(1.4,1.5))
## $root
## [1] 1.472774
## 
## $f.root
## [1] 1.661879e-05
## 
## $iter
## [1] 4
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc1LN<-uniroot(f1LN,c(1.4,1.5),tol = 1e-6)
CDFxMarketTLN(sigmaCREMR_lnt,gama1,muCREMR_lnt+log(gama1),sCREMR_lnt,xbar1,xc1LN$root)
## [1] 0.7966002
CDFxPlanTLN(sigmaCREMR_lnt,gama1,(muCREMR_lnt+log(gama1))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar1,xc1LN$root)
## [1] 0.1525987
CDFxMarketTLN(sigmaCREMR_lnt,gama1,muCREMR_lnt+log(gama1),sCREMR_lnt,xbar1,xc1LN$root)/CDFxPlanTLN(sigmaCREMR_lnt,gama1,(muCREMR_lnt+log(gama1))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar1,xc1LN$root)
## [1] 5.22023
gama3=1.5
xbar3=gama3*(omegaCREMR_lnt^sigmaCREMR_lnt+1)
f3LN=function(x){
  PDFdifferenceLN(sigmaCREMR_lnt,gama3,muCREMR_lnt+log(gama3),(muCREMR_lnt+log(gama3))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar3,x)
}
uniroot(f3LN,c(2.2,2.21))
## $root
## [1] 2.209169
## 
## $f.root
## [1] 1.706608e-05
## 
## $iter
## [1] 2
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc3LN<-uniroot(f3LN,c(2.2,2.21),tol = 1e-6)
CDFxMarketTLN(sigmaCREMR_lnt,gama3,muCREMR_lnt+log(gama3),sCREMR_lnt,xbar3,xc3LN$root)
## [1] 0.7966002
CDFxPlanTLN(sigmaCREMR_lnt,gama3,(muCREMR_lnt+log(gama3))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar3,xc3LN$root)
## [1] 0.1525987
CDFxMarketTLN(sigmaCREMR_lnt,gama3,muCREMR_lnt+log(gama3),sCREMR_lnt,xbar3,xc3LN$root)/CDFxPlanTLN(sigmaCREMR_lnt,gama3,(muCREMR_lnt+log(gama3))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar3,xc3LN$root)
## [1] 5.22023
gama5=2
xbar5=gama5*(omegaCREMR_lnt^sigmaCREMR_lnt+1)
f5LN=function(x){
  PDFdifferenceLN(sigmaCREMR_lnt,gama5,muCREMR_lnt+log(gama5),(muCREMR_lnt+log(gama5))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar5,x)
}
uniroot(f5LN,c(2.93,2.95))
## $root
## [1] 2.945531
## 
## $f.root
## [1] -2.307234e-08
## 
## $iter
## [1] 3
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc5LN<-uniroot(f5LN,c(2.93,2.95))
CDFxMarketTLN(sigmaCREMR_lnt,gama5,muCREMR_lnt+log(gama5),sCREMR_lnt,xbar5,xc5LN$root)
## [1] 0.7966002
CDFxPlanTLN(sigmaCREMR_lnt,gama5,(muCREMR_lnt+log(gama5))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar5,xc5LN$root)
## [1] 0.1525987
CDFxMarketTLN(sigmaCREMR_lnt,gama5,muCREMR_lnt+log(gama5),sCREMR_lnt,xbar5,xc5LN$root)/CDFxPlanTLN(sigmaCREMR_lnt,gama5,(muCREMR_lnt+log(gama5))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar5,xc5LN$root)
## [1] 5.22023
gama6=3
xbar6=gama6*(omegaCREMR_lnt^sigmaCREMR_lnt+1)
f6LN=function(x){
  PDFdifferenceLN(sigmaCREMR_lnt,gama6,muCREMR_lnt+log(gama6),(muCREMR_lnt+log(gama6))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar6,x)
}
uniroot(f6LN,c(4.40,4.42))
## $root
## [1] 4.418277
## 
## $f.root
## [1] -3.885796e-06
## 
## $iter
## [1] 2
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc6LN<-uniroot(f6LN,c(4.40,4.42))
CDFxMarketTLN(sigmaCREMR_lnt,gama6,muCREMR_lnt+log(gama6),sCREMR_lnt,xbar6,xc6LN$root)
## [1] 0.7966002
CDFxPlanTLN(sigmaCREMR_lnt,gama6,(muCREMR_lnt+log(gama6))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar6,xc6LN$root)
## [1] 0.1525987
CDFxMarketTLN(sigmaCREMR_lnt,gama6,muCREMR_lnt+log(gama6),sCREMR_lnt,xbar6,xc6LN$root)/CDFxPlanTLN(sigmaCREMR_lnt,gama6,(muCREMR_lnt+log(gama6))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar6,xc6LN$root)
## [1] 5.22023
gama7=4
xbar7=gama7*(omegaCREMR_lnt^sigmaCREMR_lnt+1)
f7LN=function(x){
  PDFdifferenceLN(sigmaCREMR_lnt,gama7,muCREMR_lnt+log(gama7),(muCREMR_lnt+log(gama7))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar7,x)
}
uniroot(f7LN,c(5.8,5.9))
## $root
## [1] 5.89106
## 
## $f.root
## [1] -9.493441e-08
## 
## $iter
## [1] 3
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc7LN<-uniroot(f7LN,c(5.8,5.9))
CDFxMarketTLN(sigmaCREMR_lnt,gama7,muCREMR_lnt+log(gama7),sCREMR_lnt,xbar7,xc7LN$root)
## [1] 0.7966002
CDFxPlanTLN(sigmaCREMR_lnt,gama7,(muCREMR_lnt+log(gama7))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar7,xc7LN$root)
## [1] 0.1525987
CDFxMarketTLN(sigmaCREMR_lnt,gama7,muCREMR_lnt+log(gama7),sCREMR_lnt,xbar7,xc7LN$root)/CDFxPlanTLN(sigmaCREMR_lnt,gama7,(muCREMR_lnt+log(gama7))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar7,xc7LN$root)
## [1] 5.22023
gama8=5
xbar8=gama8*(omegaCREMR_lnt^sigmaCREMR_lnt+1)
f8LN=function(x){
  PDFdifferenceLN(sigmaCREMR_lnt,gama8,muCREMR_lnt+log(gama8),(muCREMR_lnt+log(gama8))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar8,x)
}
uniroot(f8LN,c(7.3,7.4))
## $root
## [1] 7.363823
## 
## $f.root
## [1] -2.613998e-07
## 
## $iter
## [1] 3
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
xc8LN<-uniroot(f8LN,c(7.3,7.4))
CDFxMarketTLN(sigmaCREMR_lnt,gama8,muCREMR_lnt+log(gama8),sCREMR_lnt,xbar8,xc8LN$root)
## [1] 0.7966002
CDFxPlanTLN(sigmaCREMR_lnt,gama8,(muCREMR_lnt+log(gama8))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar8,xc8LN$root)
## [1] 0.1525987
CDFxMarketTLN(sigmaCREMR_lnt,gama8,muCREMR_lnt+log(gama8),sCREMR_lnt,xbar8,xc8LN$root)/CDFxPlanTLN(sigmaCREMR_lnt,gama8,(muCREMR_lnt+log(gama8))/sigmaCREMR_lnt,sCREMR_lnt,omegaCREMR_lnt,xbar8,xc8LN$root)
## [1] 5.22023