library("readxl")
my_data <- read_excel("mydatax11.xlsx")
library("readxl")
my_data <- read_excel("mydatax11.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)
# 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))
(kCREMR=estimCREMR$par[1])
(sigmaCREMR=estimCREMR$par[2])
(omegaCREMR=((sigmaCREMR-1)*mmin/(mmin+sigmaCREMR-sigmaCREMR*mmin))^(1/sigmaCREMR))
# AIC (with three estimated parameters)
pCREMR = 3
(AIC_CREMR = 2*pCREMR - 2*(-estimCREMR$value))
# MLE estimators (closed form)
(omegaLIN=2*mmin-1)
(kkLIN = nn/(sum(log(2*m-1))-nn*log(omegaLIN)))
# AIC (with two estimated parameters):
pLIN = 2
(AIC_LIN = 2*pLIN - 2*(sum(log(2*kkLIN*(omegaLIN^kkLIN)*(2*m-1)^(-kkLIN-1)))))
# MLE estimators (closed form)
(omegaLES=mmin^2)
(kkLES = 1/((2/nn)*sum(log(m))-log(omegaLES)))
# AIC (with two estimated parameters):
pLES = 2
(AIC_LES = 2*pLES - 2*(sum(log((2*kkLES*(omegaLES^(kkLES))*(m)^(-2*kkLES-1))))))
# MLE estimators (closed form)
(omegaTLOG=mmin*exp(mmin))
(kkTLOG = 1/((1/nn)*sum(log(m)+m)-log(omegaTLOG)))
# 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)))))
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)
# 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))
# Estimated parameters (untruncated lognormal)
(muCREMR_ln=estimCREMR_ln$par[1])
(zCREMR_ln=estimCREMR_ln$par[2])
(ssCREMR_ln=estimCREMR_ln$par[3])
(sigmaCREMR_ln=u(zCREMR_ln))
(sCREMR_ln=exp(ssCREMR_ln))
# Truncated lognormal case
# Truncation point
(mt=min(m))
# 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))
# Estimated parameters (truncated lognormal)
(muCREMR_lnt=estimCREMR_lnt$par[1])
(zCREMR_lnt=estimCREMR_lnt$par[2])
(ssCREMR_lnt=estimCREMR_lnt$par[3])
(sigmaCREMR_lnt=u(zCREMR_lnt))
(sCREMR_lnt=exp(ssCREMR_lnt))
(omegaCREMR_lnt=((mt*(sigmaCREMR_lnt-1))/(mt+sigmaCREMR_lnt-mt*sigmaCREMR_lnt))^(1/sigmaCREMR_lnt))
(truncCREMR_lnt=pnorm((log((mt*(sigmaCREMR_lnt-1))/(mt+sigmaCREMR_lnt-mt*sigmaCREMR_lnt))-muCREMR_lnt)/(sigmaCREMR_lnt*sCREMR_lnt),0,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-truncCREMR_lnt))
(formatC(1/(1-truncCREMR_lnt), digits = 4, format = "f"))
# AIC (with four estimated parameters):
pCREMR_lnt = 4
(AIC_CREMR_lnt = 2*pCREMR_lnt - 2*(-estimCREMR_lnt$value))
gama1=1
f1P=function(x){
PDFdifferenceP(sigmaCREMR,gama1,omegaCREMR,kCREMR, x)
}
uniroot(f1P,c(1.1,1.5))
library("readxl")
my_data <- read_excel("mydatax11.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)
# 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))
(kCREMR=estimCREMR$par[1])
(sigmaCREMR=estimCREMR$par[2])
(omegaCREMR=((sigmaCREMR-1)*mmin/(mmin+sigmaCREMR-sigmaCREMR*mmin))^(1/sigmaCREMR))
# AIC (with three estimated parameters)
pCREMR = 3
(AIC_CREMR = 2*pCREMR - 2*(-estimCREMR$value))
# MLE estimators (closed form)
(omegaLIN=2*mmin-1)
(kkLIN = nn/(sum(log(2*m-1))-nn*log(omegaLIN)))
# AIC (with two estimated parameters):
pLIN = 2
(AIC_LIN = 2*pLIN - 2*(sum(log(2*kkLIN*(omegaLIN^kkLIN)*(2*m-1)^(-kkLIN-1)))))
# MLE estimators (closed form)
(omegaLES=mmin^2)
(kkLES = 1/((2/nn)*sum(log(m))-log(omegaLES)))
# AIC (with two estimated parameters):
pLES = 2
(AIC_LES = 2*pLES - 2*(sum(log((2*kkLES*(omegaLES^(kkLES))*(m)^(-2*kkLES-1))))))
# MLE estimators (closed form)
(omegaTLOG=mmin*exp(mmin))
(kkTLOG = 1/((1/nn)*sum(log(m)+m)-log(omegaTLOG)))
# 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)))))
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)
# 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))
# Estimated parameters (untruncated lognormal)
(muCREMR_ln=estimCREMR_ln$par[1])
(zCREMR_ln=estimCREMR_ln$par[2])
(ssCREMR_ln=estimCREMR_ln$par[3])
(sigmaCREMR_ln=u(zCREMR_ln))
(sCREMR_ln=exp(ssCREMR_ln))
# Truncated lognormal case
# Truncation point
(mt=min(m))
# 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))
# Estimated parameters (truncated lognormal)
(muCREMR_lnt=estimCREMR_lnt$par[1])
(zCREMR_lnt=estimCREMR_lnt$par[2])
(ssCREMR_lnt=estimCREMR_lnt$par[3])
(sigmaCREMR_lnt=u(zCREMR_lnt))
(sCREMR_lnt=exp(ssCREMR_lnt))
(omegaCREMR_lnt=((mt*(sigmaCREMR_lnt-1))/(mt+sigmaCREMR_lnt-mt*sigmaCREMR_lnt))^(1/sigmaCREMR_lnt))
(truncCREMR_lnt=pnorm((log((mt*(sigmaCREMR_lnt-1))/(mt+sigmaCREMR_lnt-mt*sigmaCREMR_lnt))-muCREMR_lnt)/(sigmaCREMR_lnt*sCREMR_lnt),0,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-truncCREMR_lnt))
(formatC(1/(1-truncCREMR_lnt), digits = 4, format = "f"))
# AIC (with four estimated parameters):
pCREMR_lnt = 4
(AIC_CREMR_lnt = 2*pCREMR_lnt - 2*(-estimCREMR_lnt$value))
# MLE estimation (untruncated lognormal)
(muLIN=(1/nn)*sum(log(2*m-1)))
(sLIN=sqrt((1/nn)*sum((log(2*m-1)-muLIN)^2)))
# 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))
# MLE estimation (untruncated lognormal)
(muLES=(2/nn)*sum(log(m)))
(sLES=sqrt((1/nn)*sum((2*log(m)-muLES)^2)))
# 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))
# AIC (with three estimated parameters):
pLES_lnt = 3
(AIC_LES_lnt = 2*pLES_lnt - 2*(-estimLES_lnt$value))
# MLE estimation (untruncated lognormal) - closed form
(muTLOG=(1/nn)*sum(log(m*exp(m))))
(sTLOG=sqrt((1/nn)*sum((log(m*exp(m))-muTLOG)^2)))
# 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))
# Estimated parameters (truncated lognormal)
(muTLOGt=estimTLOG_lnt$par[1])
(sTLOGt=exp(estimTLOG_lnt$par[2]))
# AIC (with three estimated parameters)
pTLOG_lnt = 3
(AIC_TLOG_lnt = 2*pTLOG_lnt - 2*(-estimTLOG_lnt$value))
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))
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))
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)
}
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)
}
gama1=1
f1P=function(x){
PDFdifferenceP(sigmaCREMR,gama1,omegaCREMR,kCREMR, x)
}
uniroot(f1P,c(1.1,1.5))
xc1P<-uniroot(f1P,c(1.1,1.5),tol = 1e-6)
CDFxMarketP(sigmaCREMR,gama1,omegaCREMR,kCREMR, xc1P$root)
CDFxPlanP(sigmaCREMR,gama1,omegaCREMR,kCREMR, xc1P$root)
CDFxMarketP(sigmaCREMR,gama1,omegaCREMR,kCREMR, xc1P$root)/CDFxPlanP(sigmaCREMR,gama1,omegaCREMR,kCREMR, xc1P$root)
