import numpy as np
import pandas as pd
import numdifftools as nd
import copy
from scipy import stats
from scipy.optimize import minimize
def AppLogDen_SL79(theta, y, X, P):
# Get parameters
[alpha, beta1, beta2, beta3, sigma2u, sigma2v, = theta[:-2]
sigma2w2, sigma2w3, sigma23] = theta[-2:]
mu
= P[:, 1] - P[:, 0] + np.log(beta1) - np.log(beta2)
B2 = P[:, 2] - P[:, 0] + np.log(beta1) - np.log(beta3)
B3 = X[:, 0] - X[:, 1] - B2
w2 = X[:, 0] - X[:, 2] - B3
w3 = np.concatenate([w2.reshape(-1,1), w3.reshape(-1,1)], axis=1)
w
= np.array([[sigma2w3, sigma23],
SigmaWW
[sigma23, sigma2w2]])if np.all(np.linalg.eigvals(SigmaWW) > 0):
= beta1 + beta2 + beta3
r
= np.sqrt(sigma2u/sigma2v)
_lambda = sigma2u+sigma2v
sigma2 = np.sqrt(sigma2)
sigma
= y - alpha - X[:, 0]*beta1 - X[:, 1]*beta2 - X[:, 2]*beta3
eps = P[:, 1] - P[:, 0] + np.log(beta1) - np.log(beta2)
B2 = P[:, 2] - P[:, 0] + np.log(beta1) - np.log(beta3)
B3
= 2/sigma*stats.norm.pdf(eps/sigma, 0, 1)*stats.norm.cdf(-_lambda*eps/sigma, 0, 1)
tech_den = np.log(tech_den)
logDen_Tech
= np.log(stats.multivariate_normal.pdf(w, mean=mu, cov=SigmaWW))
logDen_Alloc
= np.log(r) + logDen_Tech + logDen_Alloc
logDen
return logDen
else:
= -1e10
logDen
return logDen
def AppLoglikelihood_SL79(theta, y, X, P):
= AppLogDen_SL79(theta, y, X, P)
logDen
= -np.sum(logDen)
logL
return logL
def AppEstimate_SL79(y, X, P):
# Starting values
= -2
alpha = 0.5
beta1 = 0.15
beta2 = 0.2;
beta3 = 0.01
sigma2u = 0.002
sigma2v = 0.01
sigma2w2 = 0.01
sigma2w3 = 0
sigma23 = 0.3
mu2 = 0.3
mu3 = [mu2, mu3]
mu
= [alpha, beta1, beta2, beta3, sigma2u, sigma2v, sigma2w2, sigma2w3,
theta0 1+sigma23)/(1-sigma23))]+ mu
np.log((
# Bounds to ensure sigma are positive
= [(None, None) for x in range(4)] + [(1e-6, np.inf) for i in range(4)] + [(None, None) for x in range(3)]
bounds
# Minimize the negative log-likelihood using numerical optimization.
= minimize(
MLE_results =AppLoglikelihood_SL79,
fun=theta0,
x0=bounds,
bounds="L-BFGS-B",
method= 1e-6,
tol ={"ftol": 1e-6, "maxiter": 1000, "maxfun": 6*1000, "maxcor": 500},
options=(y.values, X.values, P.values),
args
)
= MLE_results.x # Estimated parameter vector
theta = MLE_results.fun * -1 # Log-likelihood at the solution
logMLE
# standard errors
# gradient calculation
= 1e-6
delta = np.zeros((len(y), len(theta)))
grad for i in range(len(theta)):
= np.copy(theta)
theta1 = theta[i] + delta
theta1[i] = (AppLogDen_SL79(theta1, y.values, X.values, P.values) - AppLogDen_SL79(theta, y.values, X.values, P.values))/delta
grad[:,i]
# Hessian calculation
= nd.Hessian(f=AppLoglikelihood_SL79, method="central", step=1e-5)(theta, y.values, X.values, P.values)
hessian
= grad.T@grad
OPG = np.sqrt(np.diag(np.linalg.inv(OPG)))
ster = np.sqrt(np.diag(np.linalg.inv(-hessian)@OPG@np.linalg.inv(-hessian)))
ster_rob
return theta, ster, ster_rob, logMLE
def MatrixToCholeski(Sigma):
= np.linalg.cholesky(Sigma)
SigmaChol = SigmaChol[np.tril_indices(SigmaChol.shape[0], k=0)] # elements of lower triangle of the choleski
CholOfSigma
return CholOfSigma
def CholeskiToMatrix(CholOfSigma):
= np.zeros((3,3))
SigmaChol 0], k=0)] = CholOfSigma
SigmaChol[np.tril_indices(SigmaChol.shape[= SigmaChol@SigmaChol.T
Sigma
return Sigma
def AppLogDen_SL80(theta, y, X, P):
# Get parameters
= theta[:4]
[alpha, beta1, beta2, beta3] = theta[-3]
sigma2v = np.concatenate([np.array([0]), theta[-2:]])
mu
#Sigma = CholeskiToMatrix(theta[4:10])
= np.zeros((3,3))
Sigma 0], k=0)] = theta[4:10]
Sigma[np.tril_indices(Sigma.shape[= Sigma + Sigma.T-np.diag(np.diag(Sigma))
Sigma
= np.array(Sigma)
G 0, 0] = G[0, 0] + sigma2v
G[
= np.array(G)
Gstar 0, 1:3] = -G[0, 1:3]
Gstar[1:3, 0] = -G[1:3, 0]
Gstar[
= y - alpha - X[:,0] * beta1 - X[:,1] * beta2 - X[:,2] * beta3
eps = P[:,1] - P[:,0] + np.log(beta1) - np.log(beta2)
B2 = P[:,2] - P[:,0] + np.log(beta1) - np.log(beta3)
B3 = X[:,0] - X[:,1] - B2
w2 = X[:,0] - X[:,2] - B3
w3 = np.array([w2, w3]).T
w = np.array([eps, w2, w3]).T
all_values
= beta1 + beta2 + beta3
r
= np.linalg.inv(Sigma)
InvSigma
= (eps + sigma2v * w@InvSigma[0, 1:3]) / np.sqrt(sigma2v) * np.sqrt(np.linalg.det(Sigma) / np.linalg.det(G))
A = (eps - sigma2v * w@InvSigma[0, 1:3]) / np.sqrt(sigma2v) * np.sqrt(np.linalg.det(Sigma) / np.linalg.det(G))
Astar
= stats.multivariate_normal.pdf(all_values, mean=mu, cov=G)
term1 = stats.multivariate_normal.pdf(all_values, mean=mu, cov=Gstar)
term2 = np.log(r) + np.log(stats.norm.cdf(-A, 0, 1) * term1 + stats.norm.cdf(-Astar, 0, 1) * term2)
logDen
return logDen
def AppLoglikelihood_SL80(theta, y, X, P):
= transform_parameters(theta=theta)
theta
= AppLogDen_SL80(theta, y, X, P)
logDen
= -np.sum(logDen)
logL
return logL
def transform_parameters(theta):
= copy.deepcopy(theta)
transformed_theta = CholeskiToMatrix(theta[4:10])
Sigma 4:10] = Sigma[np.tril_indices(Sigma.shape[0], k=0)]
transformed_theta[
return transformed_theta
def AppEstimate_SL80(y, X, P):
# Starting values
= -2.6
alpha = 0.5
beta1 = 0.2
beta2 = 0.2;
beta3 = 0.14
sigma2u = 0.1
sigma2v = 0
sigmau2 = 0
sigmau3 = 0.2
sigma22 = 0.2
sigma33 = 0.1
sigma23 = 0.5
mu2 = 0.1
mu3 = [mu2, mu3]
mu
= np.array([[sigma2u, sigmau2, sigmau3],
Sigma
[sigmau2, sigma22, sigma23],
[sigmau3, sigma23, sigma33]])
= MatrixToCholeski(Sigma) # 1x6 vector of lower triangle of chol of Sigma
CholOfSigma
= np.array([alpha, beta1, beta2, beta3] + CholOfSigma.tolist() + [sigma2v] + mu)
theta0
# Bounds to ensure sigma are positive
= [(None, None) for x in range(4 + len(CholOfSigma))] + [(1e-6, np.inf) for i in range(3)]
bounds
# Minimize the negative log-likelihood using numerical optimization.
= minimize(
MLE_results =AppLoglikelihood_SL80,
fun=theta0,
x0=bounds,
bounds="Nelder-Mead",
method= 1e-6,
tol ={"maxiter": 20000},
options=(y.values, X.values, P.values),
args
)
# transform parameters
= MLE_results.x
theta_ = transform_parameters(theta=theta_)
theta = MLE_results.fun * -1 # Log-likelihood at the solution
logMLE
# standard errors
# gradient calculation
= 1e-6
delta = np.zeros((len(y), len(theta)))
grad for i in range(len(theta)):
= np.copy(theta)
theta1 = theta[i] + delta
theta1[i] = (AppLogDen_SL80(theta1, y.values, X.values, P.values) - AppLogDen_SL80(theta, y.values, X.values, P.values))/delta
grad[:,i]
# Hessian calculation
= nd.Hessian(f=AppLoglikelihood_SL80, method="central", step=1e-6)(theta_, y.values, X.values, P.values)
hessian
# Hessian transformation
= np.eye(len(theta_))
D for i in range(len(theta_)):
= copy.deepcopy(theta_)
theta1 = theta_[i] + delta
theta1[i] = (transform_parameters(theta1)-transform_parameters(theta_))/delta
D[:,i]
= grad.T@grad
OPG = np.sqrt(np.diag(np.linalg.inv(OPG)))
ster = np.sqrt(np.diag(D.T@np.linalg.inv(hessian)@D@OPG@D.T@np.linalg.inv(hessian)@D))
ster_rob
return theta, ster, ster_rob, logMLE
6 Chapter 6: Modelling Dependence in Production Systems
6.1 Exercise 6.1: The SL79 and SL80 estimators for US utilities.
This exercise implements the SL79 and SL80 estimators for the U.S. utilities dataset steamelectric.xls. The data for this exercise is the same as the data used by Rungsuriyawiboon and Stefanou (2007) and Amsler et al. (2021). The output variable is net power generated, in megawatt-hours; the inputs are fuel, labor and maintenance (L&M), and capital; the input prices are (i) the average fuel price in dollars per BTU, (ii) cost-share weighted price of labor and maintenance, where the price of labor is company-wide average wage rate and the price of maintenance is a price index of electrical supplies from the Bureau of Labor Statistics, and (iii) the yield of the latest debt issued by the firm. Quantities of fuel, L&M and capital are computed by dividing the reported fuel expenses, aggregate L&M costs and capital costs by the relevant prices. The dataset contained 1008 observations but upon data cleaning as described by Amsler et al. (2021) we ended up with 934 observations. All of the parameters are as defined earlier; in addition,
# import data
= pd.read_excel('steamelectric.xlsx')
data
'fuel_price_index'] = np.nan
data[-1, -1] = data.iloc[:-1, 2].values/data.iloc[1:,2].values
data.iloc[:'LM_price_index'] = np.nan
data[-1, -1] = data.iloc[:-1, 3].values/data.iloc[1:,3].values
data.iloc[:
13:1007:14, :] = np.nan #drop last fuel price index for each plant
data.iloc[= data.dropna()
data
= (data['Fuel Costs ($1000)']/data['fuel_price_index'])*1000*1e-6; #fuel costs over price index
X1 = data['Operating Costs ($1000)']/data['LM_price_index']*1000*1e-6; #LM costs over price index
X2 = data['Capital ($1000)']/1000; # capital
X3
= np.log(pd.concat([X1, X2, X3], axis=1))
X = np.log(data[['fuel_price_index', 'LM_price_index', 'User Cost of Capital']])
P = np.log(data['Output (MWhr)']*1e-6) # output ml MWpH
y
= AppEstimate_SL79(y=y,
(SL79_theta, SL79_sterr, SL79_ster_rob, SL79_logMLE) =X,
X=P)
P= AppEstimate_SL80(y=y,
(SL80_theta, SL80_sterr, SL80_ster_rob, SL80_logMLE) =X,
X=P)
P
# Display the results as a table
= pd.DataFrame(
results =["SL79 Est", "SL79 St.Er.(OPG)", "SL79 St.Er.(Rob)",
columns"SL80 Est", "SL80 St.Er.(OPG)", "SL80 St.Er.(Rob)"],
=[
indexr"$$\alpha$$",
r"$$\beta_{1}$$",
r"$$\beta_{2}$$",
r"$$\beta_{3}$$",
r"$$\sigma_{u}^{2}$$",
r"$$\Sigma^{2u}$$",
r"$$\Sigma^{3u}$$",
r"$$\Sigma^{22}$$",
r"$$\Sigma^{23}$$",
r"$$\Sigma^{33}$$",
r"$$\sigma^{2}_{v}$$",
r"$$\mu_{2}$$",
r"$$\mu_{3}$$",
],
)= np.round(SL79_theta.reshape(-1, 1), 3)
SL79_theta = np.round(SL79_sterr.reshape(-1, 1), 3)
SL79_sterr = np.round(SL79_ster_rob.reshape(-1, 1), 3)
SL79_ster_rob
= np.round(SL80_theta.reshape(-1, 1), 3)
SL80_theta = np.round(SL80_sterr.reshape(-1, 1), 3)
SL80_sterr = np.round(SL80_ster_rob.reshape(-1, 1), 3)
SL80_ster_rob
= np.concatenate(
SL79_frontier
[4],
SL79_theta[:4],
SL79_sterr[:4]
SL79_ster_rob[:
],=1,
axis
)= np.concatenate(
SL80_frontier
[4],
SL80_theta[:4],
SL80_sterr[:4]
SL80_ster_rob[:
],=1,
axis
)
= np.concatenate(
SL79_sigmas
[4]], np.array([[np.nan, np.nan]]).T, SL79_theta[[6]], SL79_theta[[8]], SL79_theta[[5]], SL79_theta[[9]]], axis=0),
np.concatenate([SL79_theta[[4]], np.array([[np.nan, np.nan]]).T, SL79_sterr[[6]], SL79_sterr[[8]], SL79_sterr[[5]], SL79_sterr[[9]]], axis=0),
np.concatenate([SL79_sterr[[4]], np.array([[np.nan, np.nan]]).T, SL79_ster_rob[[6]], SL79_ster_rob[[8]], SL79_ster_rob[[5]], SL79_ster_rob[[9]]], axis=0),
np.concatenate([SL79_ster_rob[[
],=1,
axis
)= np.concatenate(
SL80_sigmas
[4:6], np.array([SL80_theta[7], SL80_theta[6]]), SL80_theta[8:-2]], axis=0),
np.concatenate([SL80_theta[4:6], np.array([SL80_sterr[7], SL80_sterr[6]]), SL80_sterr[8:-2]], axis=0),
np.concatenate([SL80_sterr[4:6], np.array([SL80_ster_rob[7], SL80_ster_rob[6]]), SL80_ster_rob[8:-2]], axis=0),
np.concatenate([SL80_ster_rob[
],=1,
axis
)
= sigmas = np.concatenate(
SL79_mu
[-2:],
SL79_theta[-2:],
SL79_sterr[-2:]
SL79_ster_rob[
],=1,
axis
)= sigmas = np.concatenate(
Sl80_mu
[-2:],
SL80_theta[-2:],
SL80_sterr[-2:]
SL80_ster_rob[
],=1,
axis
)
0:4, 0:3] = SL79_frontier
results.iloc[0:4, 3:] = SL80_frontier
results.iloc[4:11, 0:3] = SL79_sigmas
results.iloc[4:11, 3:] = SL80_sigmas
results.iloc[-2:, 0:3] = SL79_mu
results.iloc[-2:, 3:] = Sl80_mu
results.iloc[ display(results)
SL79 Est | SL79 St.Er.(OPG) | SL79 St.Er.(Rob) | SL80 Est | SL80 St.Er.(OPG) | SL80 St.Er.(Rob) | |
---|---|---|---|---|---|---|
-2.613 | 0.064 | 0.073 | -2.586 | 0.064 | 0.069 | |
0.511 | 0.024 | 0.027 | 0.482 | 0.027 | 0.027 | |
0.263 | 0.013 | 0.022 | 0.267 | 0.013 | 0.021 | |
0.239 | 0.025 | 0.02 | 0.255 | 0.028 | 0.021 | |
0.167 | 0.014 | 0.014 | 0.172 | 0.014 | 0.02 | |
NaN | NaN | NaN | 0.04 | 0.015 | 0.012 | |
NaN | NaN | NaN | 0.069 | 0.012 | 0.009 | |
0.206 | 0.007 | 0.015 | 0.243 | 0.01 | 0.024 | |
0.097 | 0.006 | 0.012 | 0.093 | 0.006 | 0.017 | |
0.017 | 0.003 | 0.003 | 0.208 | 0.007 | 0.012 | |
0.568 | 0.08 | 0.123 | 0.014 | 0.003 | 0.002 | |
0.568 | 0.08 | 0.123 | 0.649 | 0.09 | 0.121 | |
-0.004 | 0.145 | 0.123 | 0.134 | 0.161 | 0.129 |
# Remove all previous function and variable definitions before the next exercise
%reset
6.2 Exercise 6.2: The APS2A and APS2B Copulas
Exercise 6.2 demonstrates how to simulate from and estimate APS2A and APS2B copulas. We generate scatterplots of 1,000 observations simulated using the APS2A and APS2B copulas, respectively. The value of
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from scipy import stats
def APS2A_deriv_u2(u1, u2, theta):
# Partial derivative of APS2A copula CDF w.r.t second argument
return u1 + theta*u1*(u1 - 1)*(4*u2**2 - 6*u2 + 2) + theta*u1*u2*(8*u2 - 6)*(u1 - 1)
def APS2A_c12_deriv_u2_inverse_func(u2, u1, theta, U):
return APS2A_deriv_u2(u2, u1, theta) - U
def APS2B_deriv_u2(u1, u2, theta):
if u2 <= 0.5:
return u1 - theta*u1*(4*u2 - 1)*(u1 - 1)
else:
return u1 + theta*u1*(4*u2 - 3)*(u1 - 1)
def APS2B_c12_deriv_u2_inverse_func(u2, u1, theta, U):
return APS2B_deriv_u2(u2, u1, theta) - U
def APS2A_loglikelihood(theta, u, w):
= stats.halfnorm.cdf(u)
F_u = stats.norm.cdf(w)
F_w
= 1 + theta*(1-2*F_u)*(1 - 12*(F_w-0.5)**2);
den = np.log(den)
logden
= np.sum(logden)*-1
LL
return LL
def APS2B_loglikelihood(theta, u, w):
= stats.halfnorm.cdf(u)
F_u = stats.norm.cdf(w)
F_w
= 1 + theta*(1-2*F_u)*(1 - 4*np.abs(F_w-0.5))
den = np.log(den)
logden
= np.sum(logden)*-1
LL
return LL
#Simulate data from APS2A Copulas
= 1000
N = 0.35
theta 432)
np.random.seed(= np.random.uniform(size=(N, 2))
U = np.zeros((N, 2))
APS2A_samples
for i in range(N):
0] = U[i,0]
APS2A_samples[i,= optimize.root_scalar(f=APS2A_c12_deriv_u2_inverse_func,
z2 =(U[i,0], theta, U[i,1]),
args=[0, 1])
bracket1] = z2.root
APS2A_samples[i,
= stats.norm.ppf(APS2A_samples[:, 0]);
APS2A_w = 1*stats.norm.ppf((APS2A_samples[:,1]+1)/2) #simulated half normal rvs
APS2A_u
0])
display(stats.pearsonr(APS2A_u, APS2A_w)[abs(APS2A_w))[0]) display(stats.pearsonr(APS2A_u, np.
0.055259235060341895
0.20849729680913273
# Simulate data from APS2B copulas
432)
np.random.seed(= 0.8
theta = np.random.uniform(size=(N, 2))
U = np.zeros((N, 2))
APS2B_samples
for i in range(N):
0] = U[i,0]
APS2B_samples[i,= optimize.root_scalar(f=APS2B_c12_deriv_u2_inverse_func,
z2 =(U[i,0], theta, U[i,1]),
args=[0, 1])
bracket1] = z2.root
APS2B_samples[i,
= stats.norm.ppf(APS2B_samples[:, 0])
APS2B_w = 1*stats.norm.ppf((APS2B_samples[:,1]+1)/2)
APS2B_u
0])
display(stats.pearsonr(APS2B_u, APS2B_w)[abs(APS2B_w))[0]) display(stats.pearsonr(APS2B_u, np.
0.05969797350051892
0.2775653676762329
# Obtain dependence parameters by fitting to data
= 0.3
theta0_APS2A
= optimize.minimize(fun=APS2A_loglikelihood,
MLE_results =theta0_APS2A,
x0="Nelder-Mead",
method=(APS2A_u, APS2A_w))
args= MLE_results.x[0]
theta_hat_APS2A = MLE_results.fun
logMLE_APS2A
= 0.7
theta0_APS2B = optimize.minimize(fun=APS2B_loglikelihood,
MLE_results =theta0_APS2B,
x0="Nelder-Mead",
method=(APS2B_u, APS2B_w))
args= MLE_results.x[0]
theta_hat_APS2B = MLE_results.fun
logMLE_APS2B
= plt.figure(figsize=(8, 8))
fig abs(APS2A_w), APS2A_u)
plt.scatter(np.r'$|\omega|$')
plt.xlabel('u')
plt.ylabel(=r'$\rho(u, \omega)$={}'.format(round(stats.pearsonr(APS2A_u, APS2A_w)[0], 3)),
plt.annotate(text=(2.2, 3.6))
xy=r'$\rho(u, |\omega|)$={}'.format(round(stats.pearsonr(APS2A_u, np.abs(APS2A_w))[0], 3)),
plt.annotate(text=(2.2, 3.4))
xy=r'$\hat\theta$={}'.format(round(theta_hat_APS2A, 3)),
plt.annotate(text=(2.2, 3.15))
xy
= plt.figure(figsize=(8, 8))
fig abs(APS2B_w), APS2B_u)
plt.scatter(np.r'$|\omega|$')
plt.xlabel('u')
plt.ylabel(=r'$\rho(u, \omega)$={}'.format(round(stats.pearsonr(APS2B_u, APS2B_w)[0], 3)),
plt.annotate(text=(2.2, 3.6))
xy=r'$\rho(u, |\omega|)$={}'.format(round(stats.pearsonr(APS2B_u, np.abs(APS2B_w))[0], 3)),
plt.annotate(text=(2.2, 3.4))
xy=r'$\hat\theta$={}'.format(round(theta_hat_APS2B, 3)),
plt.annotate(text=(2.2, 3.15)) xy
Text(2.2, 3.15, '$\\hat\\theta$=0.841')
# Remove all previous function and variable definitions before the next exercise
%reset
6.3 Exercise 6.3: The APS3A and APS3B estimators for US utilities
This exercise estimates the model in (6.1)-(6.2) using MSLE with the APS- 3-A and APS-3-B copulas on the same dataset as in Exercise 6.1. The number of replications for the simulation of the likelihood is 1,000. The results are similar to each other and to our estimates of the Schmidt and Lovell (1980) model reported in Exercise 6.1. As in Exercise 6.1, “OPG” stands for the expected outer-product of gradient and “Rob” stands for the robust version of the standard errors.
import numpy as np
import pandas as pd
import numdifftools as nd
import copy
from scipy import stats
from scipy.optimize import minimize
def MatrixToCholeski(Sigma):
= np.linalg.cholesky(Sigma)
SigmaChol = SigmaChol[np.tril_indices(SigmaChol.shape[0], k=0)] # elements of lower triangle of the choleski
CholOfSigma
return CholOfSigma
def CholeskiToMatrix(CholOfSigma):
= np.zeros((2,2))
SigmaChol 0], k=0)] = CholOfSigma
SigmaChol[np.tril_indices(SigmaChol.shape[= SigmaChol@SigmaChol.T
Sigma
return Sigma
def transform_parameters(theta):
= copy.deepcopy(theta)
transformed_theta = CholeskiToMatrix(theta[6:9])
Sigma 6:9] = Sigma[np.tril_indices(Sigma.shape[0], k=0)]
transformed_theta[
return transformed_theta
def AppLoglikelihood_APS3A_FMLE(theta, y, X, P, us):
= transform_parameters(theta=theta)
theta
= AppLogDen_APS3A_FMLE(theta, y, X, P, us)
logDen
= -np.sum(logDen)
logL
return logL
def AppLogDen_APS3A_FMLE(theta, y, X, P, us):
= len(y)
N = len(us)
S
= theta[:5]
[alpha, beta1, beta2, beta3, beta4]
= theta[4]
sigma2u = theta[5]
sigma2v
= np.zeros((2,2))
SigmaWW 0], k=0)] = theta[6:9]
SigmaWW[np.tril_indices(SigmaWW.shape[= SigmaWW + SigmaWW.T-np.diag(np.diag(SigmaWW))
SigmaWW
= theta[9]
theta12 = theta[10]
theta13 = theta[-2:]
mu
= y - alpha - X[:,0] * beta1 - X[:,1] * beta2 - X[:,2] * beta3
eps = P[:,1] - P[:,0] + np.log(beta1) - np.log(beta2)
B2 = P[:,2] - P[:,0] + np.log(beta1) - np.log(beta3)
B3 = X[:,0] - X[:,1] - B2
w2 = X[:,0] - X[:,2] - B3
w3 = np.array([w2, w3]).T
w
= stats.norm.pdf(w2, mu[0], np.sqrt(SigmaWW[0,0]))
DenW2 = stats.norm.pdf(w3, mu[1], np.sqrt(SigmaWW[1,1]))
DenW3 = stats.norm.cdf(w2, mu[0], np.sqrt(SigmaWW[0,0]))
CdfW2 = stats.norm.cdf(w3, mu[1], np.sqrt(SigmaWW[1,1]))
CdfW3
= 2*(stats.norm.cdf(np.sqrt(sigma2u)*us, 0, np.sqrt(sigma2u)) -0.5);
CdfUs
= np.repeat(eps.reshape(-1,1), S, axis=1).T
eps_SxN = np.repeat(np.sqrt(sigma2u)*us.reshape(-1,1), N, axis=1)
us_SxN
= eps_SxN + us_SxN
EpsPlusUs
= stats.norm.pdf(EpsPlusUs, 0, np.sqrt(sigma2v));
DenEpsPlusUs
= np.repeat(CdfW2.reshape(-1,1), S, axis=1).T
w2s_SxN = np.repeat(CdfW3.reshape(-1,1), S, axis=1).T
w3s_SxN = np.repeat(CdfUs, N, axis=1)
w1s_SxN = 1 + theta12*(1-2*w1s_SxN)*(1-12*(w2s_SxN -0.5)**2)
c12_SxN = 1 + theta13*(1-2*w1s_SxN)*(1-12*(w3s_SxN -0.5)**2)
c13_SxN = stats.multivariate_normal.pdf(w, mean=mu, cov=SigmaWW)/(DenW2*DenW3)
c23 = np.repeat(c23.reshape(-1,1), S, axis=1).T
c23_SxN
= 1+(c12_SxN-1)+(c13_SxN-1)+(c23_SxN-1)
cstar
= np.mean(cstar*DenEpsPlusUs, axis=0)
Integral = DenW2*DenW3*Integral
DenAll
= beta1+beta2+beta3;
r
= np.log(r)+np.log(DenAll)
logDen
return logDen
def AppEstimate_APS3A_FMLE(y, X, P):
= 1000 # No. simulated draws for the integral
S 12345)
np.random.seed(= stats.norm.ppf((np.random.uniform(size=(S, 1))+1)/2, 0, 1)
us
= -2.6
alpha = 0.5
beta1 = 0.2
beta2 = 0.2
beta3 = 0.14
sigma2u = 0.01
sigma2v = -0.1
theta12 = -0.1
theta13 = 0.337
sigma22 = 0.59
sigma33 = 0.20
sigma23 = 0.3
mu2 = 0.3
mu3 = [mu2, mu3]
mu
= np.array([[sigma22, sigma23],
SigmaWW
[sigma23, sigma33]])= MatrixToCholeski(SigmaWW) # vector of lower triange of chol of SigmaWW
CholOfSigma
= np.array([alpha, beta1, beta2, beta3, sigma2u, sigma2v] + CholOfSigma.tolist() + [theta12, theta13] + mu)
theta0
# Bounds to ensure sigma are positive
= [(None, None) for x in range(4)] + [(1e-6, np.inf) for i in range(2)] + [(None, None) for x in range(len(CholOfSigma)+4)]
bounds
# Minimize the negative log-likelihood using numerical optimization.
= minimize(AppLoglikelihood_APS3A_FMLE,
MLE_results
theta0,="Nelder-Mead",
method= 1e-6,
tol =bounds,
bounds={"maxiter": 3000},
options=(y.values, X.values, P.values, us))
args
# transform parameters
= MLE_results.x
theta_ = transform_parameters(theta=theta_)
theta = MLE_results.fun * -1 # Log-likelihood at the solution
logMLE
# standard errors
# gradient calculation
= 1e-6
delta = np.zeros((len(y), len(theta)))
grad for i in range(len(theta)):
= np.copy(theta)
theta1 = theta[i] + delta
theta1[i] = (AppLogDen_APS3A_FMLE(theta1, y.values, X.values, P.values, us) - AppLogDen_APS3A_FMLE(theta, y.values, X.values, P.values, us))/delta
grad[:,i]
# Hessian calculation
= nd.Hessian(f=AppLoglikelihood_APS3A_FMLE, method="central", step=1e-6)(theta_, y.values, X.values, P.values, us)
hessian
# Hessian transformation
= np.eye(len(theta_))
D for i in range(len(theta_)):
= copy.deepcopy(theta_)
theta1 = theta_[i] + delta
theta1[i] = (transform_parameters(theta1)-transform_parameters(theta_))/delta
D[:,i] = grad.T@grad
OPG
# Delta method
= np.sqrt(np.diag(np.linalg.inv(OPG)))
ster_OPG = np.sqrt(np.diag(D.T@np.linalg.inv(hessian)@D@OPG@D.T@np.linalg.inv(hessian)@D))
ster_rob
return theta, ster_OPG, ster_rob, logMLE
def AppLoglikelihood_APS3B_FMLE(theta, y, X, P, us):
= transform_parameters(theta=theta)
theta
= AppLogDen_APS3B_FMLE(theta, y, X, P, us)
logDen
= -np.sum(logDen)
logL
return logL
def AppLogDen_APS3B_FMLE(theta, y, X, P, us):
= len(y)
N = len(us)
S
= theta[:5]
[alpha, beta1, beta2, beta3, beta4]
= theta[4]
sigma2u = theta[5]
sigma2v
= np.zeros((2,2))
SigmaWW 0], k=0)] = theta[6:9]
SigmaWW[np.tril_indices(SigmaWW.shape[= SigmaWW + SigmaWW.T-np.diag(np.diag(SigmaWW))
SigmaWW
= theta[9]
theta12 = theta[10]
theta13 = theta[-2:]
mu
= y - alpha - X[:,0] * beta1 - X[:,1] * beta2 - X[:,2] * beta3
eps = P[:,1] - P[:,0] + np.log(beta1) - np.log(beta2)
B2 = P[:,2] - P[:,0] + np.log(beta1) - np.log(beta3)
B3 = X[:,0] - X[:,1] - B2
w2 = X[:,0] - X[:,2] - B3
w3 = np.array([w2, w3]).T
w
= stats.norm.pdf(w2, mu[0], np.sqrt(SigmaWW[0,0]))
DenW2 = stats.norm.pdf(w3, mu[1], np.sqrt(SigmaWW[1,1]))
DenW3 = stats.norm.cdf(w2, mu[0], np.sqrt(SigmaWW[0,0]))
CdfW2 = stats.norm.cdf(w3, mu[1], np.sqrt(SigmaWW[1,1]))
CdfW3
= 2*(stats.norm.cdf(np.sqrt(sigma2u)*us, 0, np.sqrt(sigma2u)) -0.5);
CdfUs
= np.repeat(eps.reshape(-1,1), S, axis=1).T
eps_SxN = np.repeat(np.sqrt(sigma2u)*us.reshape(-1,1), N, axis=1)
us_SxN
= eps_SxN + us_SxN
EpsPlusUs
= stats.norm.pdf(EpsPlusUs, 0, np.sqrt(sigma2v));
DenEpsPlusUs
= np.repeat(CdfW2.reshape(-1,1), S, axis=1).T
w2s_SxN = np.repeat(CdfW3.reshape(-1,1), S, axis=1).T
w3s_SxN = np.repeat(CdfUs, N, axis=1)
w1s_SxN = 1 + theta12*(1-2*w1s_SxN)*(1-4*np.abs(w2s_SxN -0.5))
c12_SxN = 1 + theta13*(1-2*w1s_SxN)*(1-4*np.abs(w3s_SxN -0.5))
c13_SxN = stats.multivariate_normal.pdf(w, mean=mu, cov=SigmaWW)/(DenW2*DenW3)
c23 = np.repeat(c23.reshape(-1,1), S, axis=1).T
c23_SxN
= 1+(c12_SxN-1)+(c13_SxN-1)+(c23_SxN-1)
cstar
= np.mean(cstar*DenEpsPlusUs, axis=0)
Integral = DenW2*DenW3*Integral
DenAll
= beta1+beta2+beta3;
r
= np.log(r)+np.log(DenAll)
logDen
return logDen
def AppEstimate_APS3B_FMLE(y, X, P):
= 1000 # No. simulated draws for the integral
S 12345)
np.random.seed(= stats.norm.ppf((np.random.uniform(size=(S, 1))+1)/2, 0, 1)
us
= -2.6
alpha = 0.5
beta1 = 0.2
beta2 = 0.2
beta3 = 0.14
sigma2u = 0.01
sigma2v = -0.1
theta12 = -0.1
theta13 = 0.337
sigma22 = 0.59
sigma33 = 0.20
sigma23 = 0.3
mu2 = 0.3
mu3 = [mu2, mu3]
mu
= np.array([[sigma22, sigma23],
SigmaWW
[sigma23, sigma33]])= MatrixToCholeski(SigmaWW) # vector of lower triange of chol of SigmaWW
CholOfSigma
= np.array([alpha, beta1, beta2, beta3, sigma2u, sigma2v] + CholOfSigma.tolist() + [theta12, theta13] + mu)
theta0
# Bounds to ensure sigma are positive
= [(None, None) for x in range(4)] + [(1e-6, np.inf) for i in range(2)] + [(None, None) for x in range(len(CholOfSigma)+4)]
bounds
# Minimize the negative log-likelihood using numerical optimization.
= minimize(AppLoglikelihood_APS3B_FMLE,
MLE_results
theta0,="Nelder-Mead",
method= 1e-6,
tol =bounds,
bounds={"maxiter": 3000},
options=(y.values, X.values, P.values, us))
args
# transform parameters
= MLE_results.x
theta_ = transform_parameters(theta=theta_)
theta = MLE_results.fun * -1 # Log-likelihood at the solution
logMLE
# standard errors
# gradient calculation
= 1e-6
delta = np.zeros((len(y), len(theta)))
grad for i in range(len(theta)):
= np.copy(theta)
theta1 = theta[i] + delta
theta1[i] = (AppLogDen_APS3B_FMLE(theta1, y.values, X.values, P.values, us) - AppLogDen_APS3B_FMLE(theta, y.values, X.values, P.values, us))/delta
grad[:,i]
# Hessian calculation
= nd.Hessian(f=AppLoglikelihood_APS3B_FMLE, method="central", step=1e-6)(theta_, y.values, X.values, P.values, us)
hessian
# Hessian transformation
= np.eye(len(theta_))
D for i in range(len(theta_)):
= copy.deepcopy(theta_)
theta1 = theta_[i] + delta
theta1[i] = (transform_parameters(theta1)-transform_parameters(theta_))/delta
D[:,i] = grad.T@grad
OPG
# Delta method
= np.sqrt(np.diag(np.linalg.inv(OPG)))
ster_OPG = np.sqrt(np.diag(D.T@np.linalg.inv(hessian)@D@OPG@D.T@np.linalg.inv(hessian)@D))
ster_rob
return theta, ster_OPG, ster_rob, logMLE
# import data
= pd.read_excel('steamelectric.xlsx')
data
'fuel_price_index'] = np.nan
data[-1, -1] = data.iloc[:-1, 2].values/data.iloc[1:,2].values
data.iloc[:'LM_price_index'] = np.nan
data[-1, -1] = data.iloc[:-1, 3].values/data.iloc[1:,3].values
data.iloc[:
13:1007:14, :] = np.nan #drop last fuel price index for each plant
data.iloc[= data.dropna()
data
= (data['Fuel Costs ($1000)']/data['fuel_price_index'])*1000*1e-6; #fuel costs over price index
X1 = data['Operating Costs ($1000)']/data['LM_price_index']*1000*1e-6; #LM costs over price index
X2 = data['Capital ($1000)']/1000; # capital
X3
= np.log(pd.concat([X1, X2, X3], axis=1))
X = np.log(data[['fuel_price_index', 'LM_price_index', 'User Cost of Capital']])
P = np.log(data['Output (MWhr)']*1e-6) # output ml MWpH y
# import data
= pd.read_excel('steamelectric.xlsx')
data
'fuel_price_index'] = np.nan
data[-1, -1] = data.iloc[:-1, 2].values/data.iloc[1:,2].values
data.iloc[:'LM_price_index'] = np.nan
data[-1, -1] = data.iloc[:-1, 3].values/data.iloc[1:,3].values
data.iloc[:
13:1007:14, :] = np.nan #drop last fuel price index for each plant
data.iloc[= data.dropna()
data
= (data['Fuel Costs ($1000)']/data['fuel_price_index'])*1000*1e-6; #fuel costs over price index
X1 = data['Operating Costs ($1000)']/data['LM_price_index']*1000*1e-6; #LM costs over price index
X2 = data['Capital ($1000)']/1000; # capital
X3
= np.log(pd.concat([X1, X2, X3], axis=1))
X = np.log(data[['fuel_price_index', 'LM_price_index', 'User Cost of Capital']])
P = np.log(data['Output (MWhr)']*1e-6) # output ml MWpH
y
= AppEstimate_APS3A_FMLE(y=y, X=X, P=P)
APS3A_theta, APS3A_ster_OPG, APS3A_ster_rob, APS3A_logMLE = AppEstimate_APS3B_FMLE(y=y, X=X, P=P)
APS3B_theta, APS3B_ster_OPG, APS3B_ster_rob, APS3B_logMLE
# Display the results as a table
= pd.DataFrame(
results =["AS3A Est", "AS3A St.Er.(OPG)", "AS3A St.Er.(Rob)",
columns"AS3B Est", "AS3B St.Er.(OPG)", "AS3B St.Er.(Rob)"],
=[
indexr"$$\alpha$$",
r"$$\beta_{1}$$",
r"$$\beta_{2}$$",
r"$$\beta_{3}$$",
r"$$\sigma_{u}^{2}$$",
r"$$\Sigma^{22}$$",
r"$$\Sigma^{23}$$",
r"$$\Sigma^{33}$$",
r"$$\sigma^{2}_{v}$$",
r"$$\mu_{2}$$",
r"$$\mu_{3}$$",
r"$$\theta_{12}$$",
r"$$\theta_{13}$$"
],
)= np.round(APS3A_theta.reshape(-1, 1), 3)
APS3A_theta = np.round(APS3A_ster_OPG.reshape(-1, 1), 3)
APS3A_ster_OPG = np.round(APS3A_ster_rob.reshape(-1, 1), 3)
APS3A_ster_rob
= np.round(APS3B_theta.reshape(-1, 1), 3)
APS3B_theta = np.round(APS3B_ster_OPG.reshape(-1, 1), 3)
APS3B_ster_OPG = np.round(APS3B_ster_rob.reshape(-1, 1), 3)
APS3B_ster_rob
= np.concatenate(
APS3A_frontier
[4],
APS3A_theta[:4],
APS3A_ster_OPG[:4]
APS3A_ster_rob[:
],=1,
axis
)= np.concatenate(
APS3B_frontier
[4],
APS3B_theta[:4],
APS3B_ster_OPG[:4]
APS3B_ster_rob[:
],=1,
axis
)
= np.concatenate(
APS3A_sigmas
[4]], APS3A_theta[6:9], APS3A_theta[[5]]], axis=0),
np.concatenate([APS3A_theta[[4]], APS3A_ster_OPG[6:9], APS3A_ster_OPG[[5]]], axis=0),
np.concatenate([APS3A_ster_OPG[[4]], APS3A_ster_rob[6:9], APS3A_ster_rob[[5]]], axis=0),
np.concatenate([APS3A_ster_rob[[
],=1,
axis
)
= np.concatenate(
APS3B_sigmas
[4]], APS3B_theta[6:9], APS3B_theta[[5]]], axis=0),
np.concatenate([APS3B_theta[[4]], APS3B_ster_OPG[6:9], APS3B_ster_OPG[[5]]], axis=0),
np.concatenate([APS3B_ster_OPG[[4]], APS3B_ster_rob[6:9], APS3B_ster_rob[[5]]], axis=0),
np.concatenate([APS3B_ster_rob[[
],=1,
axis
)
= np.concatenate(
APS3A_thetas
[-4:-2],
APS3A_theta[-4:-2],
APS3A_ster_OPG[-4:-2]
APS3A_ster_rob[
],=1,
axis
)
= np.concatenate(
APS3B_thetas
[-4:-2],
APS3B_theta[-4:-2],
APS3B_ster_OPG[-4:-2]
APS3B_ster_rob[
],=1,
axis
)
= sigmas = np.concatenate(
APS3A_mu
[-2:],
APS3A_theta[-2:],
APS3A_ster_OPG[-2:]
APS3A_ster_rob[
],=1,
axis
)= sigmas = np.concatenate(
APS3B_mu
[-2:],
APS3B_theta[-2:],
APS3B_ster_OPG[-2:]
APS3B_ster_rob[
],=1,
axis
)
0:4, 0:3] = APS3A_frontier
results.iloc[0:4, 3:] = APS3B_frontier
results.iloc[4:9, 0:3] = APS3A_sigmas
results.iloc[4:9, 3:] = APS3B_sigmas
results.iloc[9:11, 0:3] = APS3A_mu
results.iloc[9:11, 3:] = APS3B_mu
results.iloc[-2:, 0:3] = APS3A_thetas
results.iloc[-2:, 3:] = APS3B_thetas
results.iloc[ display(results)
/var/folders/q0/9kb46g7n1_57dfk7vtdsv8t4w169ch/T/ipykernel_8109/1901131725.py:96: RuntimeWarning: invalid value encountered in log
logDen = np.log(r)+np.log(DenAll)
/var/folders/q0/9kb46g7n1_57dfk7vtdsv8t4w169ch/T/ipykernel_8109/1901131725.py:235: RuntimeWarning: invalid value encountered in log
logDen = np.log(r)+np.log(DenAll)
AS3A Est | AS3A St.Er.(OPG) | AS3A St.Er.(Rob) | AS3B Est | AS3B St.Er.(OPG) | AS3B St.Er.(Rob) | |
---|---|---|---|---|---|---|
-2.733 | 0.068 | 0.1 | -2.731 | 0.062 | 0.069 | |
0.544 | 0.025 | 0.039 | 0.493 | 0.025 | 0.028 | |
0.218 | 0.016 | 0.041 | 0.253 | 0.014 | 0.022 | |
0.256 | 0.025 | 0.017 | 0.279 | 0.024 | 0.018 | |
0.151 | 0.016 | 0.019 | 0.18 | 0.015 | 0.017 | |
0.249 | 0.011 | 0.021 | 0.247 | 0.01 | 0.02 | |
0.115 | 0.007 | 0.012 | 0.119 | 0.007 | 0.012 | |
0.219 | 0.008 | 0.013 | 0.209 | 0.007 | 0.009 | |
0.024 | 0.005 | 0.006 | 0.018 | 0.004 | 0.003 | |
0.306 | 0.099 | 0.251 | 0.562 | 0.089 | 0.132 | |
-0.001 | 0.134 | 0.1 | 0.181 | 0.128 | 0.109 | |
-0.198 | 0.032 | 0.034 | -0.28 | 0.062 | 0.057 | |
-0.035 | 0.058 | 0.065 | -0.088 | 0.077 | 0.077 |
# Remove all previous function and variable definitions before the next exercise
%reset
6.4 Exercise 6.4 Vine APS-2-copula estimates for the 111 US utilities.
This exercise estimates the model in (6.1)-(6.2) using MSLE with the vine copulas based on the APS-2-A and APS-2-B copulas. The data file is cowing.xlsx which was also used in Chapter 2. In this exercise we use an approach known as the extended value extension. That is, rather than use a parameter transformation of the correlation matrix we check if the correlation matrix is positive definite within the objective function and return a very large value of the objective if this condition is not met.
Table 6.3 reports the results. We also report the estimates based on the APS3 copula for this dataset. The copula dependence parameter estimates are, for the most part, large in magnitude but statistically insignificant for this example.
import numpy as np
import pandas as pd
import copy
from scipy import stats
from scipy.optimize import minimize
def AppLoglikelihood_APS3A_FMLE(theta, y, X, P, us):
= AppLogDen_APS3A_FMLE(theta, y, X, P, us)
logDen
= -np.sum(logDen)
logL
return logL
def AppLogDen_APS3A_FMLE(theta, y, X, P, us):
= len(y)
n = len(us)
S
[alpha, beta1, beta2, beta3, mu2, mu3, sigma2u, sigma2v, = theta[:]
sigma2w2, sigma2w3, pw1w2, theta12, theta13]
= y - alpha - X[:,0] * beta1 - X[:,1] * beta2 - X[:,2] * beta3
eps = P[:,1] - P[:,0] + np.log(beta1) - np.log(beta2)
B2 = P[:,2] - P[:,0] + np.log(beta1) - np.log(beta3)
B3 = X[:,0] - X[:,1] - B2
w2 = X[:,0] - X[:,2] - B3
w3
= 2*(stats.norm.cdf(np.sqrt(sigma2u)*us, 0, np.sqrt(sigma2u)) -0.5)
CdfUs = np.repeat(eps.reshape(-1,1), S, axis=1).T
eps_SxN
= np.sqrt(sigma2u)*us
us_SxN = eps_SxN + us_SxN
EpsPlusUs = stats.norm.pdf(EpsPlusUs, 0, np.sqrt(sigma2v))
DenEpsPlusUs
=np.array([[1,pw1w2],
RhoWW 1]])
[pw1w2,
# Check Rho
if not np.all(np.linalg.eigvals(RhoWW) > 0):
= np.ones((n, 1))*-1e8
logDen else:
# Evaluate the integral via simulation (to integrate out u from eps+u)
= np.zeros(us_SxN.shape)
simulated_copula_pdfs = stats.norm.pdf(w2, mu2, np.sqrt(sigma2w2))
DenW2 = stats.norm.pdf(w3, mu3, np.sqrt(sigma2w3))
DenW3
= np.repeat(w2.reshape(-1,1), S, axis=1).T
w2_rep = np.repeat(w3.reshape(-1,1), S, axis=1).T
w3_rep # Compute the CDF (standard normal) for repeated allocative inefficiency terms
= stats.norm.cdf(w2_rep, mu2, np.sqrt(sigma2w2))
CDF_w2_rep = stats.norm.cdf(w3_rep, mu3, np.sqrt(sigma2w3))
CDF_w3_rep
for j in range(n):
= 1 + theta12*(1-2*CdfUs[:,j])*(1 - 12*(CDF_w2_rep[:,j] - 0.5)**2)
simulated_c12_pdf = 1 + theta13*(1-2*CdfUs[:,j])*(1 - 12*(CDF_w3_rep[:,j] - 0.5)**2)
simulated_c13_pdf = np.concatenate([stats.norm.ppf(CDF_w2_rep[:, j]).reshape(-1,1),
U -1,1)],
stats.norm.ppf(CDF_w3_rep[:, j]).reshape(=1)
axis= stats.multivariate_normal.pdf(U, mean = np.array([0, 0]), cov=RhoWW)/ np.prod(stats.norm.pdf(U), axis=1)
simulated_c23_pdf = 1 + (simulated_c12_pdf -1) + (simulated_c13_pdf - 1) + (simulated_c23_pdf - 1)
APS3A_copula_pdf = APS3A_copula_pdf
simulated_copula_pdfs[:,j]
= np.mean(simulated_copula_pdfs*DenEpsPlusUs, axis=0)
Integral
= DenW2*DenW3*Integral
DenAll
= beta1 + beta2 + beta3
r = np.log(r) + np.log(DenAll)
logDen
return logDen
def AppEstimate_APS3A_FMLE(y, X, P):
= len(y)
n = 100 # No. simulated draws for the integral
S 12345)
np.random.seed(= stats.norm.ppf((np.random.uniform(size=(S, n))+1)/2, 0, 1)
us
# Starting values for MLE
= -11.6
initial_alpha = 0.06
initial_beta1 = 1.04
initial_beta2 = 0.06
initial_beta3 = 0.014
initial_sigma2u = 0.002
initial_sigma2v = 0 # APS c12 copula parameter
initial_theta12 = 0 # APS c13 copula parameter
initial_theta13 = 0.33 # variance of w2
initial_sigma2w2 = 0.59 # variance of w3
initial_sigma2w3 = 1.3 # mean of normal distribution for w1
initial_mu2 = 1 #mean of normal distribution for w2
initial_mu3
# Find starting value for rho in gaussian copula for w1w2
= P.iloc[:,1] - P.iloc[:,0] + np.log(initial_beta1) - np.log(initial_beta2)
B2 = P.iloc[:,2] - P.iloc[:,0] + np.log(initial_beta1) - np.log(initial_beta3)
B3 = X.iloc[:, 0] - X.iloc[:, 1] - B2
w2 = X.iloc[:, 0] - X.iloc[:, 2] - B3
w3
= stats.kendalltau(x=stats.norm.cdf(w2, initial_mu2, np.sqrt(initial_sigma2w2)),
initial_pw1w2 =stats.norm.cdf(w3, initial_mu3, np.sqrt(initial_sigma2w2)))[0]
y
= [initial_alpha, initial_beta1, initial_beta2, initial_beta3,
theta0
initial_mu2, initial_mu3,
initial_sigma2u, initial_sigma2v, initial_sigma2w2,
initial_sigma2w3, initial_pw1w2, initial_theta12, initial_theta13]
# Bounds to ensure sigma2v and sigma2u are >= 0
= [(None, None) for x in range(6)] + [
bounds 1e-5, np.inf),
(1e-5, np.inf),
(+ [(None, None) for x in range(5)]
]
# Minimize the negative log-likelihood using numerical optimization.
= minimize(
MLE_results =AppLoglikelihood_APS3A_FMLE,
fun=theta0,
x0="L-BFGS-B",
method= 1e-6,
tol ={"ftol": 1e-6, "maxiter": 1000, "maxfun": 6*1000},
options=(y.values, X.values, P.values, us),
args=bounds,
bounds
)
= MLE_results.x
theta = MLE_results.fun * -1 # Log-likelihood at the solution
logMLE
# Estimation of standard errors
= 1e-6;
delta = np.zeros((len(y), len(theta)))
grad for i in range(len(theta)):
= copy.deepcopy(theta)
theta1 = theta[i] + delta
theta1[i] = (AppLogDen_APS3A_FMLE(theta1, y.values, X.values, P.values, us) - AppLogDen_APS3A_FMLE(theta, y.values, X.values, P.values, us))/delta
grad[:,i]
= grad.T@grad # Outer product of the gradient
OPG = np.sqrt(np.diag(np.linalg.inv(OPG)))
sterr
return theta, sterr, logMLE
def AppLoglikelihood_APS3B_FMLE(theta, y, X, P, us):
= AppLogDen_APS3B_FMLE(theta, y, X, P, us)
logDen
= -np.sum(logDen)
logL
return logL
def AppLogDen_APS3B_FMLE(theta, y, X, P, us):
= len(y)
n = len(us)
S
[alpha, beta1, beta2, beta3, mu2, mu3, sigma2u, sigma2v, = theta[:]
sigma2w2, sigma2w3, pw1w2, theta12, theta13]
= y - alpha - X[:,0] * beta1 - X[:,1] * beta2 - X[:,2] * beta3
eps = P[:,1] - P[:,0] + np.log(beta1) - np.log(beta2)
B2 = P[:,2] - P[:,0] + np.log(beta1) - np.log(beta3)
B3 = X[:,0] - X[:,1] - B2
w2 = X[:,0] - X[:,2] - B3
w3
= 2*(stats.norm.cdf(np.sqrt(sigma2u)*us, 0, np.sqrt(sigma2u)) -0.5)
CdfUs = np.repeat(eps.reshape(-1,1), S, axis=1).T
eps_SxN
= np.sqrt(sigma2u)*us
us_SxN = eps_SxN + us_SxN
EpsPlusUs = stats.norm.pdf(EpsPlusUs, 0, np.sqrt(sigma2v))
DenEpsPlusUs
=np.array([[1,pw1w2],
RhoWW 1]])
[pw1w2,
# Check Rho
if not np.all(np.linalg.eigvals(RhoWW) > 0):
= np.ones((n, 1))*-1e8
logDen else:
# Evaluate the integral via simulation (to integrate out u from eps+u)
= np.zeros(us_SxN.shape)
simulated_copula_pdfs = stats.norm.pdf(w2, mu2, np.sqrt(sigma2w2))
DenW2 = stats.norm.pdf(w3, mu3, np.sqrt(sigma2w3))
DenW3
= np.repeat(w2.reshape(-1,1), S, axis=1).T
w2_rep = np.repeat(w3.reshape(-1,1), S, axis=1).T
w3_rep # Compute the CDF (standard normal) for repeated allocative inefficiency terms
= stats.norm.cdf(w2_rep, mu2, np.sqrt(sigma2w2))
CDF_w2_rep = stats.norm.cdf(w3_rep, mu3, np.sqrt(sigma2w3))
CDF_w3_rep
for j in range(n):
= 1 + theta12*(1-2*CdfUs[:,j])*(1 - 4*abs(CDF_w2_rep[:,j]-0.5))
simulated_c12_pdf = 1 + theta13*(1-2*CdfUs[:,j])*(1-4*abs(CDF_w3_rep[:,j] - 0.5))
simulated_c13_pdf = np.concatenate([stats.norm.ppf(CDF_w2_rep[:, j]).reshape(-1,1),
U -1,1)],
stats.norm.ppf(CDF_w3_rep[:, j]).reshape(=1)
axis= stats.multivariate_normal.pdf(U, mean = np.array([0, 0]), cov=RhoWW)/ np.prod(stats.norm.pdf(U), axis=1)
simulated_c23_pdf = 1 + (simulated_c12_pdf -1) + (simulated_c13_pdf - 1) + (simulated_c23_pdf - 1)
APS3A_copula_pdf = APS3A_copula_pdf
simulated_copula_pdfs[:,j]
= np.mean(simulated_copula_pdfs*DenEpsPlusUs, axis=0)
Integral
= DenW2*DenW3*Integral
DenAll
= beta1 + beta2 + beta3
r = np.log(r) + np.log(DenAll)
logDen
return logDen
def AppEstimate_APS3B_FMLE(y, X, P):
= len(y)
n = 100 # No. simulated draws for the integral
S 12345)
np.random.seed(= stats.norm.ppf((np.random.uniform(size=(S, n))+1)/2, 0, 1)
us
# Starting values for MLE
= -11.6
initial_alpha = 0.06
initial_beta1 = 1.04
initial_beta2 = 0.06
initial_beta3 = 0.014
initial_sigma2u = 0.002
initial_sigma2v = 0 # APS c12 copula parameter
initial_theta12 = 0 # APS c13 copula parameter
initial_theta13 = 0.33 # variance of w2
initial_sigma2w2 = 0.59 # variance of w3
initial_sigma2w3 = 1.3 # mean of normal distribution for w1
initial_mu2 = 1 #mean of normal distribution for w2
initial_mu3
# Find starting value for rho in gaussian copula for w1w2
= P.iloc[:,1] - P.iloc[:,0] + np.log(initial_beta1) - np.log(initial_beta2)
B2 = P.iloc[:,2] - P.iloc[:,0] + np.log(initial_beta1) - np.log(initial_beta3)
B3 = X.iloc[:, 0] - X.iloc[:, 1] - B2
w2 = X.iloc[:, 0] - X.iloc[:, 2] - B3
w3
= stats.kendalltau(x=stats.norm.cdf(w2, initial_mu2, np.sqrt(initial_sigma2w2)),
initial_pw1w2 =stats.norm.cdf(w3, initial_mu3, np.sqrt(initial_sigma2w2)))[0]
y
= [initial_alpha, initial_beta1, initial_beta2, initial_beta3,
theta0
initial_mu2, initial_mu3,
initial_sigma2u, initial_sigma2v, initial_sigma2w2,
initial_sigma2w3, initial_pw1w2, initial_theta12, initial_theta13]
# Bounds to ensure sigma2v and sigma2u are >= 0
= [(None, None) for x in range(6)] + [
bounds 1e-5, np.inf),
(1e-5, np.inf),
(+ [(None, None) for x in range(5)]
]
# Minimize the negative log-likelihood using numerical optimization.
= minimize(
MLE_results =AppLoglikelihood_APS3B_FMLE,
fun=theta0,
x0="L-BFGS-B",
method= 1e-6,
tol ={"ftol": 1e-6, "maxiter": 1000, "maxfun": 6*1000},
options=(y.values, X.values, P.values, us),
args=bounds,
bounds
)
= MLE_results.x
theta = MLE_results.fun * -1 # Log-likelihood at the solution
logMLE
# Estimation of standard errors
= 1e-6;
delta = np.zeros((len(y), len(theta)))
grad for i in range(len(theta)):
= copy.deepcopy(theta)
theta1 = theta[i] + delta
theta1[i] = (AppLogDen_APS3B_FMLE(theta1, y.values, X.values, P.values, us) - AppLogDen_APS3B_FMLE(theta, y.values, X.values, P.values, us))/delta
grad[:,i]
= grad.T@grad # Outer product of the gradient
OPG = np.sqrt(np.diag(np.linalg.inv(OPG)))
sterr
return theta, sterr, logMLE
def AppLogDen_Vine_APS2A(theta, y, X, P, us):
= len(y)
N = us.shape
S, p
[alpha, beta1, beta2, beta3, mu2, mu3, sigma2u, sigma2v, = theta[:]
sigma2w2, sigma2w3, pw1w2, theta12, theta13]
= y - alpha - X[:, 0]*beta1 - X[:, 1]*beta2 - X[:, 2]*beta3
eps = P[:,1] - P[:,0] + np.log(beta1) - np.log(beta2)
B2 = P[:,2] - P[:,0] + np.log(beta1) - np.log(beta3)
B3 = X[:,0] - X[:,1] - B2
w2 = X[:,0] - X[:,2] - B3
w3 = 2*(stats.norm.cdf(np.sqrt(sigma2u)*us, 0, np.sqrt(sigma2u)) -0.5)
CdfUs = np.repeat(eps.reshape(-1,1), S, axis=1).T
eps_SxN
= np.sqrt(sigma2u)*us
us_SxN = eps_SxN + us_SxN
EpsPlusUs = stats.norm.pdf(EpsPlusUs, 0, np.sqrt(sigma2v))
DenEpsPlusUs
= stats.norm.pdf(w2, mu2, np.sqrt(sigma2w2))
DenW2 = stats.norm.pdf(w3, mu3, np.sqrt(sigma2w3))
DenW3 =np.array([[1,pw1w2],
RhoWW 1]])
[pw1w2,
# Check Rho
if not np.all(np.linalg.eigvals(RhoWW) > 0):
= np.ones((N, 1))*-1e8
logDen else:
# Evaluate the integral via simulation (to integrate out u from eps+u)
= np.zeros(us_SxN.shape)
simulated_copula_pdfs = np.repeat(w2.reshape(-1,1), S, axis=1).T
w2_rep = np.repeat(w3.reshape(-1,1), S, axis=1).T
w3_rep # Compute the CDF (standard normal) for repeated allocative inefficiency terms
= stats.norm.cdf(w2_rep, mu2, np.sqrt(sigma2w2))
CDF_w2_rep = stats.norm.cdf(w3_rep, mu3, np.sqrt(sigma2w3))
CDF_w3_rep
= CDF_w2_rep + theta12*CdfUs*CDF_w2_rep*(4*CDF_w2_rep**2 - 6*CDF_w2_rep + 2) + theta12*CDF_w2_rep*(CdfUs - 1)*(4*CDF_w2_rep**2 - 6*CDF_w2_rep + 2)
conditional_CDF_w2_u = CDF_w3_rep + theta13*CdfUs*CDF_w3_rep*(4*CDF_w3_rep**2 - 6*CDF_w3_rep + 2) + theta13*CDF_w3_rep*(CdfUs - 1)*(4*CDF_w3_rep**2 - 6*CDF_w3_rep + 2)
conditional_CDF_w3_u for j in range(N):
= 1 + theta12*(1-2*CdfUs[:,j])*(1 - 12*(CDF_w2_rep[:,j] - 0.5)**2)
simulated_c21_pdf = 1 + theta13*(1-2*CdfUs[:,j])*(1 - 12*(CDF_w3_rep[:,j] - 0.5)**2)
simulated_c31_pdf = np.concatenate([stats.norm.ppf(conditional_CDF_w2_u[:, j]).reshape(-1,1),
U -1,1)],
stats.norm.ppf(conditional_CDF_w3_u[:, j]).reshape(=1)
axis= stats.multivariate_normal.pdf(U, mean = np.array([0, 0]), cov=RhoWW)/ np.prod(stats.norm.pdf(U), axis=1)
simulated_c23_pdf = simulated_c21_pdf*simulated_c31_pdf*simulated_c23_pdf
APS2A_copula_pdf = APS2A_copula_pdf
simulated_copula_pdfs[:,j] = np.mean(simulated_copula_pdfs*DenEpsPlusUs, axis=0)
Integral
= DenW2*DenW3*Integral
DenAll
= beta1 + beta2 + beta3
r = np.log(r) + np.log(DenAll)
logDen
return logDen
def AppLoglikelihood_vine_APS2A(theta, y, X, P, us):
= AppLogDen_Vine_APS2A(theta, y, X, P, us)
logDen
= -np.sum(logDen)
logL
return logL
def AppEstimate_Vine_APS2A(y, X, P):
# Standard Normal random variables for simulation of the integral.
12345)
np.random.seed(= 100
S = len(y)
n = stats.norm.ppf((np.random.uniform(size=(S, n))+1)/2, 0, 1)
us
# Starting values for MLE
= -11.6
initial_alpha = 0.06
initial_beta1 = 1.04
initial_beta2 = 0.06
initial_beta3 = 0.014
initial_sigma2u = 0.002
initial_sigma2v = 0.2 # APS c12 copula parameter
initial_theta12 = -0.2 # APS c13 copula parameter
initial_theta13 = 0.33 # variance of w2
initial_sigma2w2 = 0.59 # variance of w3
initial_sigma2w3 = 2 # mean of normal distribution for w1
initial_mu2 = 0.5 #mean of normal distribution for w2
initial_mu3
# Find starting value for rho in gaussian copula for w1w2
= P.iloc[:,1] - P.iloc[:,0] + np.log(initial_beta1) - np.log(initial_beta2)
B2 = P.iloc[:,2] - P.iloc[:,0] + np.log(initial_beta1) - np.log(initial_beta3)
B3 = X.iloc[:, 0] - X.iloc[:, 1] - B2
w2 = X.iloc[:, 0] - X.iloc[:, 2] - B3
w3
# Conditional CDF transforms for allocative inefficiency correlation starting value
= 2*(stats.norm.cdf(np.sqrt(initial_sigma2u)*us[0, :], 0, np.sqrt(initial_sigma2u)) -0.5)
CDF_simulated_us = stats.norm.cdf(w2, initial_mu2, np.sqrt(initial_sigma2w2))
CDF_w2 = stats.norm.cdf(w3, initial_mu3 , np.sqrt(initial_sigma2w3))
CDF_w3 = CDF_w2 + (initial_theta12*CDF_w2)*(1 - CDF_w2)*(2*(6*CDF_simulated_us - 6*CDF_simulated_us**2 - 1))
APS2A_CDF_w2_conditional_u = CDF_w3 + (initial_theta13*CDF_w3)*(1 - CDF_w3)*(2*(6*CDF_simulated_us - 6*CDF_simulated_us**2 - 1))
APS2A_CDF_w3_conditional_u = stats.kendalltau(x=APS2A_CDF_w2_conditional_u,
initial_pw1w2 =APS2A_CDF_w3_conditional_u)[0]
y
= [initial_alpha, initial_beta1, initial_beta2, initial_beta3,
theta0
initial_mu2, initial_mu3, initial_sigma2u, initial_sigma2v,
initial_sigma2w2, initial_sigma2w3, initial_pw1w2,
initial_theta12, initial_theta13]
# Bounds to ensure sigma2v and sigma2u are >= 0
= [(None, None) for x in range(6)] + [
bounds 1e-5, np.inf),
(1e-5, np.inf),
(+ [(None, None) for x in range(5)]
]
# Minimize the negative log-likelihood using numerical optimization.
= minimize(
MLE_results =AppLoglikelihood_vine_APS2A,
fun=theta0,
x0="L-BFGS-B",
method= 1e-6,
tol ={"ftol": 1e-6, "maxiter": 1000, "maxfun": 6*1000},
options=(y.values, X.values, P.values, us),
args=bounds,
bounds
)
= MLE_results.x # Estimated parameter vector
theta = MLE_results.fun * -1 # Log-likelihood at the solution
log_likelihood
# Estimation of standard errors
= 1e-6;
delta = np.zeros((len(y), len(theta)))
grad for i in range(len(theta)):
= copy.deepcopy(theta)
theta1 = theta[i] + delta
theta1[i] = (AppLogDen_Vine_APS2A(theta1, y.values, X.values, P.values, us) - AppLogDen_Vine_APS2A(theta, y.values, X.values, P.values, us))/delta
grad[:,i]
= grad.T@grad # Outer product of the gradient
OPG = np.sqrt(np.diag(np.linalg.inv(OPG)))
sterr
return theta, sterr, log_likelihood
def AppLogDen_Vine_APS2B(theta, y, X, P, us):
= len(y)
n = us.shape
S, p
[alpha, beta1, beta2, beta3, mu2, mu3, sigma2u, sigma2v, = theta[:]
sigma2w2, sigma2w3, pw1w2, theta12, theta13]
= y - alpha - X[:, 0]*beta1 - X[:, 1]*beta2 - X[:, 2]*beta3
eps = P[:,1] - P[:,0] + np.log(beta1) - np.log(beta2)
B2 = P[:,2] - P[:,0] + np.log(beta1) - np.log(beta3)
B3 = X[:,0] - X[:,1] - B2
w2 = X[:,0] - X[:,2] - B3
w3 = 2*(stats.norm.cdf(np.sqrt(sigma2u)*us, 0, np.sqrt(sigma2u)) -0.5)
CdfUs = np.repeat(eps.reshape(-1,1), S, axis=1).T
eps_SxN
= np.sqrt(sigma2u)*us
us_SxN = eps_SxN + us_SxN
EpsPlusUs = stats.norm.pdf(EpsPlusUs, 0, np.sqrt(sigma2v))
DenEpsPlusUs
= stats.norm.pdf(w2, mu2, np.sqrt(sigma2w2))
DenW2 = stats.norm.pdf(w3, mu3, np.sqrt(sigma2w3))
DenW3 =np.array([[1,pw1w2],
RhoWW 1]])
[pw1w2,
# Check Rho
if not np.all(np.linalg.eigvals(RhoWW) > 0):
= np.ones((n, 1))*-1e8
logDen else:
# Evaluate the integral via simulation (to integrate out u from eps+u)
= np.zeros(us_SxN.shape)
simulated_copula_pdfs = np.repeat(w2.reshape(-1,1), S, axis=1).T
w2_rep = np.repeat(w3.reshape(-1,1), S, axis=1).T
w3_rep # Compute the CDF (standard normal) for repeated allocative inefficiency terms
= stats.norm.cdf(w2_rep, mu2, np.sqrt(sigma2w2))
CDF_w2_rep = stats.norm.cdf(w3_rep, mu3, np.sqrt(sigma2w3))
CDF_w3_rep #Compute the conditional CDFs for allocative inefficiency terms F(w1|u), F(w2|u)
= np.zeros((S, n))
conditional_CDF_w2_u = np.zeros((S, n))
conditional_CDF_w3_u
for j in range(n):
if w2[j] <= 0.5:
= CDF_w2_rep[:,j] + theta12*CdfUs[:,j]*(-2*CDF_w2_rep[:,j]**2 + CDF_w2_rep[:,j]) + theta12*(-2*CDF_w2_rep[:,j]**2 + CDF_w2_rep[:,j])*(CdfUs[:,j] - 1)
conditional_CDF_w2_u[:,j] else:
= CDF_w2_rep[:,j] + theta12*(CdfUs[:,j] - 1)*(2*CDF_w2_rep[:,j]**2 - 3*CDF_w2_rep[:,j] + 1) + theta12*CdfUs[:,j]*(2*CDF_w2_rep[:,j]**2 - 3*CDF_w2_rep[:,j] + 1)
conditional_CDF_w2_u[:,j] if w3[j] <= 0.5:
= CDF_w3_rep[:,j] + theta13*CdfUs[:,j]*(-2*CDF_w3_rep[:,j]**2 + CDF_w3_rep[:,j]) + theta13*(-2*CDF_w3_rep[:,j]**2 + CDF_w3_rep[:,j])*(CdfUs[:,j] - 1)
conditional_CDF_w3_u[:,j] else:
= CDF_w3_rep[:,j] + theta13*(CdfUs[:,j] - 1)*(2*CDF_w3_rep[:,j]**2 - 3*CDF_w3_rep[:,j] + 1) + theta13*CdfUs[:,j]*(2*CDF_w3_rep[:,j]**2 - 3*CDF_w3_rep[:,j] + 1)
conditional_CDF_w3_u[:,j]
= 1 + theta12*(1-2*CdfUs[:,j])*(1 - 4*np.abs(CDF_w2_rep[:,j] - 0.5))
simulated_c21_pdf = 1 + theta13*(1-2*CdfUs[:,j])*(1-4*np.abs(CDF_w3_rep[:,j] - 0.5))
simulated_c31_pdf = np.concatenate([stats.norm.ppf(conditional_CDF_w2_u[:, j]).reshape(-1,1),
U -1,1)],
stats.norm.ppf(conditional_CDF_w3_u[:, j]).reshape(=1)
axis= stats.multivariate_normal.pdf(U, mean = np.array([0, 0]), cov=RhoWW)/ np.prod(stats.norm.pdf(U), axis=1)
simulated_c23_pdf = simulated_c21_pdf*simulated_c31_pdf*simulated_c23_pdf
APS2B_copula_pdf
= APS2B_copula_pdf
simulated_copula_pdfs[:,j] = np.mean(simulated_copula_pdfs*DenEpsPlusUs, axis=0)
Integral
= DenW2*DenW3*Integral
DenAll
= beta1 + beta2 + beta3
r = np.log(r) + np.log(DenAll)
logDen
return logDen
def AppLoglikelihood_vine_APS2B(theta, y, X, P, us):
= AppLogDen_Vine_APS2B(theta, y, X, P, us)
logDen
= -np.sum(logDen)
logL
return logL
def AppEstimate_Vine_APS2B(y, X, P):
# Standard Normal random variables for simulation of the integral.
12345)
np.random.seed(= 100
S = len(y)
n = stats.norm.ppf((np.random.uniform(size=(S, n))+1)/2, 0, 1)
us
# Starting values for MLE
= -11.6
initial_alpha = 0.045
initial_beta1 = 1.04
initial_beta2 = 0.06
initial_beta3 = 0.014
initial_sigma2u = 0.002
initial_sigma2v = 0.1 # APS c12 copula parameter
initial_theta12 = -0.2 # APS c13 copula parameter
initial_theta13 = 0.33 # variance of w2
initial_sigma2w2 = 0.59 # variance of w3
initial_sigma2w3 = 2 # mean of normal distribution for w1
initial_mu2 = 0.5 #mean of normal distribution for w2
initial_mu3
# Find starting value for rho in gaussian copula for w1w2
= P.iloc[:,1] - P.iloc[:,0] + np.log(initial_beta1) - np.log(initial_beta2)
B2 = P.iloc[:,2] - P.iloc[:,0] + np.log(initial_beta1) - np.log(initial_beta3)
B3 = X.iloc[:, 0] - X.iloc[:, 1] - B2
w2 = X.iloc[:, 0] - X.iloc[:, 2] - B3
w3
# Conditional CDF transforms for allocative inefficiency correlation starting value
= 2*(stats.norm.cdf(np.sqrt(initial_sigma2u)*us[0, :], 0, np.sqrt(initial_sigma2u)) -0.5)
CDF_simulated_us = stats.norm.cdf(w2, initial_mu2, np.sqrt(initial_sigma2w2))
CDF_w2 = stats.norm.cdf(w3, initial_mu3 , np.sqrt(initial_sigma2w3))
CDF_w3
= CDF_w2
APS2B_CDF_w2_conditional_u = np.where(APS2B_CDF_w2_conditional_u <= 0.5)[0]
idx1 = CDF_w2 + (initial_theta12*CDF_w2)*(1 - CDF_w2)*(4*CDF_simulated_us - 1)
term1 = np.where(APS2B_CDF_w2_conditional_u > 0.5)[0]
idx2 = CDF_w2 + (initial_theta12*CDF_w2)*(1 - CDF_w2)*(3 - 4*CDF_simulated_us)
term2 = term1[idx1]
APS2B_CDF_w2_conditional_u[idx1] = term2[idx2]
APS2B_CDF_w2_conditional_u[idx2]
= CDF_w3;
APS2B_CDF_w3_conditional_u = np.where(APS2B_CDF_w3_conditional_u <= 0.5)
idx1 = CDF_w3 + (initial_theta13*CDF_w3)*(1 - CDF_w3)*(4*CDF_simulated_us - 1)
term1 = np.where(APS2B_CDF_w3_conditional_u > 0.5)
idx2 = CDF_w3 + (initial_theta12*CDF_w3)*(1 - CDF_w3)*(3 - 4*CDF_simulated_us)
term2 = term1[idx1]
APS2B_CDF_w3_conditional_u[idx1] = term2[idx2]
APS2B_CDF_w3_conditional_u[idx2]
= stats.kendalltau(x=APS2B_CDF_w2_conditional_u,
initial_pw1w2 =APS2B_CDF_w3_conditional_u)[0]
y
= [initial_alpha, initial_beta1, initial_beta2, initial_beta3,
theta0
initial_mu2, initial_mu3, initial_sigma2u, initial_sigma2v,
initial_sigma2w2, initial_sigma2w3, initial_pw1w2,
initial_theta12, initial_theta13]
# Bounds to ensure sigma2v and sigma2u are >= 0
= [(None, None) for x in range(6)] + [
bounds 1e-5, np.inf),
(1e-5, np.inf),
(+ [(None, None) for x in range(5)]
]
# Minimize the negative log-likelihood using numerical optimization.
= minimize(
MLE_results =AppLoglikelihood_vine_APS2B,
fun=theta0,
x0="Nelder-Mead",
method= 1e-6,
tol ={"maxiter": 2000},
options=(y.values, X.values, P.values, us),
args=bounds,
bounds
)
= MLE_results.x # Estimated parameter vector
theta = MLE_results.fun * -1 # Log-likelihood at the solution
log_likelihood
# Estimation of standard errors
= 1e-6;
delta = np.zeros((len(y), len(theta)))
grad for i in range(len(theta)):
= copy.deepcopy(theta)
theta1 = theta[i] + delta
theta1[i] = (AppLogDen_Vine_APS2B(theta1, y.values, X.values, P.values, us) - AppLogDen_Vine_APS2B(theta, y.values, X.values, P.values, us))/delta
grad[:,i]
= grad.T@grad # Outer product of the gradient
OPG = np.sqrt(np.diag(np.linalg.inv(OPG)))
sterr
return theta, sterr, log_likelihood
# import data
= pd.read_excel('cowing.xlsx') # Inputs and outputs are in Logarithms
production_data
= production_data['y'] # Output
y = production_data[['X1', 'X2', 'X3']]
X = np.log(production_data[['P1', 'P2', 'P3']])
P
= AppEstimate_APS3A_FMLE(y=y, X=X, P=P)
APS3A_theta, APS3A_ster, APS3A_logMLE = AppEstimate_APS3B_FMLE(y=y, X=X, P=P)
APS3B_theta, APS3B_ster, APS3B_logMLE = AppEstimate_Vine_APS2A(y=y, X=X, P=P)
vine_APS2A_theta, vine_APS2A_ster, vine_APS2A_logMLE = AppEstimate_Vine_APS2B(y=y, X=X, P=P)
vine_APS2B_theta, vine_APS2B_ster, vine_APS2B_logMLE
# Display the results as a table
= pd.DataFrame(
results =["AS3A Est", "AS3A St.Err.",
columns"AS3B Est", "AS3B St.Err",
"Vine AS2A Est", "Vine AS2A St.Err.",
"Vine AS2B Est", "Vine AS2B St.Err"],
=[
indexr"$$\alpha$$",
r"$$\alpha_{1}$$",
r"$$\alpha_{2}$$",
r"$$\alpha_{3}$$",
r"$$\mu_{1}$$",
r"$$\mu_{2}$$",
r"$$\sigma_{u}^{2}$$",
r"$$\sigma^{2}_{v}$$",
r"$$\Sigma^{2}_{1}$$",
r"$$\Sigma^{2}_{2}$$",
r"$$\rho_{12}$$",
r"$$\theta_{12}$$",
r"$$\theta_{13}$$"
],
)= np.round(APS3A_theta.reshape(-1, 1), 3)
APS3A_theta = np.round(APS3A_ster.reshape(-1, 1), 3)
APS3A_ster = np.round(APS3B_theta.reshape(-1, 1), 3)
APS3B_theta = np.round(APS3B_ster.reshape(-1, 1), 3)
APS3B_ster = np.round(vine_APS2A_theta.reshape(-1, 1), 3)
vine_APS2A_theta = np.round(vine_APS2A_ster.reshape(-1, 1), 3)
vine_APS2A_ster = np.round(vine_APS2B_theta.reshape(-1, 1), 3)
vine_APS2B_theta = np.round(vine_APS2B_ster.reshape(-1, 1), 3)
vine_APS2B_ster
= np.concatenate(
APS_frontier
[4], APS3A_ster[:4],
APS3A_theta[:4], APS3B_ster[:4],
APS3B_theta[:4], vine_APS2A_ster[:4],
vine_APS2A_theta[:4], vine_APS2B_ster[:4],
vine_APS2B_theta[:
],=1,
axis
)
= np.concatenate(
APS_sigmas
[6:10], APS3A_ster[6:10],
APS3A_theta[6:10], APS3B_ster[6:10],
APS3B_theta[6:10], vine_APS2A_ster[6:10],
vine_APS2A_theta[6:10], vine_APS2B_ster[6:10],
vine_APS2B_theta[
],=1,
axis
)
= np.concatenate(
APS_mu
[4:6], APS3A_ster[4:6],
APS3A_theta[4:6], APS3B_ster[4:6],
APS3B_theta[4:6], vine_APS2A_ster[4:6],
vine_APS2A_theta[4:6], vine_APS2B_ster[4:6],
vine_APS2B_theta[
],=1,
axis
)
= np.concatenate(
APS_thetas
[-3:], APS3A_ster[-3:],
APS3A_theta[-3:], APS3B_ster[-3:],
APS3B_theta[-3:], vine_APS2A_ster[-3:],
vine_APS2A_theta[-3:], vine_APS2B_ster[-3:],
vine_APS2B_theta[
],=1,
axis
)
0:4, :] = APS_frontier
results.iloc[4:6, :] = APS_mu
results.iloc[6:10, :] = APS_sigmas
results.iloc[-3:, :] = APS_thetas
results.iloc[
display(results)print('APS3A LL: ', APS3A_logMLE)
print('APS3B LL: ', APS3B_logMLE)
print('Vine APS2A LL: ', vine_APS2A_logMLE)
print('Vine APS2B LL: ', vine_APS2B_logMLE)
/var/folders/q0/9kb46g7n1_57dfk7vtdsv8t4w169ch/T/ipykernel_8109/1741302720.py:70: RuntimeWarning: invalid value encountered in log
logDen = np.log(r) + np.log(DenAll)
AS3A Est | AS3A St.Err. | AS3B Est | AS3B St.Err | Vine AS2A Est | Vine AS2A St.Err. | Vine AS2B Est | Vine AS2B St.Err | |
---|---|---|---|---|---|---|---|---|
-11.471 | 0.241 | -11.381 | 0.236 | -11.502 | 0.268 | -11.453 | 0.262 | |
0.045 | 0.025 | 0.041 | 0.024 | 0.045 | 0.026 | 0.058 | 0.022 | |
1.073 | 0.025 | 1.076 | 0.025 | 1.079 | 0.026 | 1.059 | 0.022 | |
0.03 | 0.029 | 0.025 | 0.029 | 0.026 | 0.03 | 0.031 | 0.029 | |
1.951 | 0.591 | 2.037 | 0.64 | 1.96 | 0.61 | 1.668 | 0.415 | |
0.662 | 1.17 | 0.588 | 1.374 | 0.508 | 1.327 | 0.444 | 0.985 | |
0.01 | 0.003 | 0.01 | 0.003 | 0.01 | 0.004 | 0.015 | 0.003 | |
0.004 | 0.001 | 0.003 | 0.001 | 0.003 | 0.001 | 0.002 | 0.0 | |
0.333 | 0.032 | 0.337 | 0.034 | 0.343 | 0.035 | 0.349 | 0.036 | |
0.582 | 0.094 | 0.604 | 0.103 | 0.593 | 0.099 | 0.575 | 0.093 | |
0.508 | 0.1 | 0.533 | 0.097 | 0.483 | 0.101 | 0.412 | 0.107 | |
0.644 | 0.368 | 0.639 | 0.458 | 0.206 | 0.415 | 0.106 | 0.453 | |
-0.17 | 0.329 | -0.089 | 0.459 | -0.22 | 0.349 | -0.118 | 0.411 |
APS3A LL: -72.52373198800693
APS3B LL: -73.12293492860093
Vine APS2A LL: -74.4857292736222
Vine APS2B LL: -77.29226482218182
# Remove all previous function and variable definitions before the next exercise
%reset