/* Mata function library for
Ryo Okui and Takahide Yanagi. Kernel Estimation for Panel Data
with Heterogeneous Dynamics. 2019.
Ryo Okui and Takahide Yanagi. Panel Data Analysis with Heterogeneous
Dynamics. 2019. */
/* Contents
1. Basic
1.1. Autocovariance
1.2. Autocorrelation
2. Empirical CDF Estimation
2.1. Naive Estimation
2.2. Half-Panel-Jackknife Estimation
2.3. Third-Order-Jackknife Estimation
3. Moment Estimation
3.1. Naive Estimation
3.2. Half-Panel-Jackknife Estimation
3.3. Third-Order-Jackknife Estimation
4. Kernel Density Estimation
4.1. Naive Estimation
4.2. Half-Panel-Jackknife Estimation
4.3. Third-Order-Jackknife Estimation
*/
version 14.0
mata:
mata clear
/// 1. Basic
// 1.1. Autocovariance
function mataacov (y, acov_order){
k = acov_order
s = length(y)
mean_est = sum(y) / s
y1 = y[1::s-k]
y2 = y[k+1::s]
acov_est = sum((y1 :- mean_est) :* (y2 :- mean_est)) / (s - k)
return(acov_est)
}
// 1.2. Autocorrelation
function mataacor (y, acor_order){
acor_est = mataacov(y, acor_order) / mataacov(y, 0)
return(acor_est)
}
/// 2. Empirical CDF estimation
function m_freq(x, X){
est = mean(X:<=x)
return(est)
}
// 2.1. Naive Estimation
function m_neecdf(data, acov_order, acor_order, B) {
N = rows(data)
S = cols(data)
grid = 100
mean_est = J(N, 1, 0)
acov_est = J(N, 1, 0)
acor_est = J(N, 1, 0)
for (i = 1; i <= N; i++) {
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i=1; i<=grid; i++) {
mean_ecdf[i] = m_freq(mean_grid[i], mean_est)
acov_ecdf[i] = m_freq(acov_grid[i], acov_est)
acor_ecdf[i] = m_freq(acor_grid[i], acor_est)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = m_freq(mean_grid[i], mean_est[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = m_freq(acov_grid[i], acov_est[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = m_freq(acor_grid[i], acor_est[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
temp=st_addvar("double", "mean_ecdf")
temp=st_addvar("double", "acov_ecdf")
temp=st_addvar("double", "acor_ecdf")
temp=st_addvar("double", "mean_grid")
temp=st_addvar("double", "acov_grid")
temp=st_addvar("double", "acor_grid")
temp=st_addvar("double", "mean_LCI")
temp=st_addvar("double", "acov_LCI")
temp=st_addvar("double", "acor_LCI")
temp=st_addvar("double", "mean_UCI")
temp=st_addvar("double", "acov_UCI")
temp=st_addvar("double", "acor_UCI")
st_addobs(max((0,grid - st_nobs())))
st_store(.,"mean_ecdf", mean_ecdf\J(st_nobs()-rows(mean_ecdf),1,.))
st_store(.,"acov_ecdf", acov_ecdf\J(st_nobs()-rows(acov_ecdf),1,.))
st_store(.,"acor_ecdf", acor_ecdf\J(st_nobs()-rows(acor_ecdf),1,.))
st_store(.,"mean_grid", mean_grid\J(st_nobs()-rows(mean_grid),1,.))
st_store(.,"acov_grid", acov_grid\J(st_nobs()-rows(acov_grid),1,.))
st_store(.,"acor_grid", acor_grid\J(st_nobs()-rows(acor_grid),1,.))
st_store(.,"mean_LCI", mean_LCI\J(st_nobs()-rows(mean_LCI),1,.))
st_store(.,"acov_LCI", acov_LCI\J(st_nobs()-rows(acov_LCI),1,.))
st_store(.,"acor_LCI", acor_LCI\J(st_nobs()-rows(acor_LCI),1,.))
st_store(.,"mean_UCI", mean_UCI\J(st_nobs()-rows(mean_UCI),1,.))
st_store(.,"acov_UCI", acov_UCI\J(st_nobs()-rows(acov_UCI),1,.))
st_store(.,"acor_UCI", acor_UCI\J(st_nobs()-rows(acor_UCI),1,.))
}
// 2.2. Half-Panel-Jackknife Estimation
function hpjecdfest1(x, X, X1, X2) {
est = m_freq(x, X)
est1 = m_freq(x, X1)
est2 = m_freq(x, X2)
hpjest = 2 * est - (est1 + est2) / 2
if (hpjest < 0) {
hpjest = 0
}
if (hpjest > 1){
hpjest = 1
}
return(hpjest)
}
function hpjecdfest2(x, X, X1, X2, X3, X4) {
est = m_freq(x, X)
est1 = m_freq(x, X1)
est2 = m_freq(x, X2)
est3 = m_freq(x, X3)
est4 = m_freq(x, X4)
hpjest = 2 * est - (est1 + est2 + est3 + est4) / 4
if (hpjest < 0) {
hpjest = 0
}
if (hpjest > 1){
hpjest = 1
}
return(hpjest)
}
function m_hpjecdf(data, acov_order, acor_order, B) {
N = rows(data)
S = cols(data)
grid = 100
if (mod(S,2)==0) {
mean_est = J(N, 1, 0)
acov_est = J(N, 1, 0)
acor_est = J(N, 1, 0)
data1 = data[,1::(S/2)]
data2 = data[,(S/2+1)::S]
mean_est1 = J(N, 1, 0)
mean_est2 = J(N, 1, 0)
acov_est1 = J(N, 1, 0)
acor_est1 = J(N, 1, 0)
acov_est2 = J(N, 1, 0)
acor_est2 = J(N, 1, 0)
for(i=1; i<=N ;i++){
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
mean_est1[i] = mean(data1[i,.]')
acov_est1[i] = mataacov(data1[i,.], acov_order)
acor_est1[i] = mataacor(data1[i,.], acor_order)
mean_est2[i] = mean(data2[i,.]')
acov_est2[i] = mataacov(data2[i,.], acov_order)
acor_est2[i] = mataacor(data2[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i=1; i<= grid; i++){
mean_ecdf[i] = hpjecdfest1(mean_grid[i], mean_est, mean_est1, mean_est2)
acov_ecdf[i] = hpjecdfest1(acov_grid[i], acov_est, acov_est1, acov_est2)
acor_ecdf[i] = hpjecdfest1(acor_grid[i], acor_est, acor_est1, acor_est2)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = hpjecdfest1(mean_grid[i], mean_est[index_boot], mean_est1[index_boot], mean_est2[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = hpjecdfest1(acov_grid[i], acov_est[index_boot], acov_est1[index_boot], acov_est2[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = hpjecdfest1(acor_grid[i], acor_est[index_boot], acor_est1[index_boot], acor_est2[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
mean_data = sort((mean_ecdf, mean_LCI, mean_UCI), 1)
acov_data = sort((acov_ecdf, acov_LCI, acov_UCI), 1)
acor_data = sort((acor_ecdf, acor_LCI, acor_UCI), 1)
mean_ecdf = mean_data[,1]
mean_LCI = mean_data[,2]
mean_UCI = mean_data[,3]
acov_ecdf = acov_data[,1]
acov_LCI = acov_data[,2]
acov_UCI = acov_data[,3]
acor_ecdf = acor_data[,1]
acor_LCI = acor_data[,2]
acor_UCI = acor_data[,3]
}
else{
mean_est = J(N, 1, 0)
acov_est = J(N, 1, 0)
acor_est = J(N, 1, 0)
data1 = data[., 1::floor(S/2)]
data2 = data[., (floor(S/2)+1)::S]
data3 = data[.,1::(floor(S/2)+1)]
data4 = data[.,(floor(S/2)+2)::S]
mean_est1 = J(N, 1, 0)
mean_est2 = J(N, 1, 0)
mean_est3 = J(N, 1, 0)
mean_est4 = J(N, 1, 0)
acov_est1 = J(N, 1, 0)
acov_est2 = J(N, 1, 0)
acov_est3 = J(N, 1, 0)
acov_est4 = J(N, 1, 0)
acor_est1 = J(N, 1, 0)
acor_est2 = J(N, 1, 0)
acor_est3 = J(N, 1, 0)
acor_est4 = J(N, 1, 0)
for(i=1; i<=N; i++){
mean_est[i] = mean(data[i,.]')
mean_est1[i] = mean(data1[i,.]')
mean_est2[i] = mean(data2[i,.]')
mean_est3[i] = mean(data3[i,.]')
mean_est4[i] = mean(data4[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acov_est1[i] = mataacov(data1[i,.], acov_order)
acov_est2[i] = mataacov(data2[i,.], acov_order)
acov_est3[i] = mataacov(data3[i,.], acov_order)
acov_est4[i] = mataacov(data4[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
acor_est1[i] = mataacor(data1[i,.], acor_order)
acor_est2[i] = mataacor(data2[i,.], acor_order)
acor_est3[i] = mataacor(data3[i,.], acor_order)
acor_est4[i] = mataacor(data4[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i=1; i<= grid; i++){
mean_ecdf[i] = hpjecdfest2(mean_grid[i], mean_est, mean_est1, mean_est2, mean_est3, mean_est4)
acov_ecdf[i] = hpjecdfest2(acov_grid[i], acov_est, acov_est1, acov_est2, acov_est3, acov_est4)
acor_ecdf[i] = hpjecdfest2(acor_grid[i], acor_est, acor_est1, acor_est2, acor_est3, acor_est4)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = hpjecdfest2(mean_grid[i], mean_est[index_boot], mean_est1[index_boot], mean_est2[index_boot], mean_est3[index_boot], mean_est4[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = hpjecdfest2(acov_grid[i], acov_est[index_boot], acov_est1[index_boot], acov_est2[index_boot], acov_est3[index_boot], acov_est4[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = hpjecdfest2(acor_grid[i], acor_est[index_boot], acor_est1[index_boot], acor_est2[index_boot], acor_est3[index_boot], acor_est4[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
mean_data = sort((mean_ecdf, mean_LCI, mean_UCI), 1)
acov_data = sort((acov_ecdf, acov_LCI, acov_UCI), 1)
acor_data = sort((acor_ecdf, acor_LCI, acor_UCI), 1)
mean_ecdf = mean_data[,1]
mean_LCI = mean_data[,2]
mean_UCI = mean_data[,3]
acov_ecdf = acov_data[,1]
acov_LCI = acov_data[,2]
acov_UCI = acov_data[,3]
acor_ecdf = acor_data[,1]
acor_LCI = acor_data[,2]
acor_UCI = acor_data[,3]
}
temp=st_addvar("double", "mean_ecdf")
temp=st_addvar("double", "acov_ecdf")
temp=st_addvar("double", "acor_ecdf")
temp=st_addvar("double", "mean_grid")
temp=st_addvar("double", "acov_grid")
temp=st_addvar("double", "acor_grid")
temp=st_addvar("double", "mean_LCI")
temp=st_addvar("double", "acov_LCI")
temp=st_addvar("double", "acor_LCI")
temp=st_addvar("double", "mean_UCI")
temp=st_addvar("double", "acov_UCI")
temp=st_addvar("double", "acor_UCI")
st_addobs(max((0,grid - st_nobs())))
st_store(.,"mean_ecdf", mean_ecdf\J(st_nobs()-rows(mean_ecdf),1,.))
st_store(.,"acov_ecdf", acov_ecdf\J(st_nobs()-rows(acov_ecdf),1,.))
st_store(.,"acor_ecdf", acor_ecdf\J(st_nobs()-rows(acor_ecdf),1,.))
st_store(.,"mean_grid", mean_grid\J(st_nobs()-rows(mean_grid),1,.))
st_store(.,"acov_grid", acov_grid\J(st_nobs()-rows(acov_grid),1,.))
st_store(.,"acor_grid", acor_grid\J(st_nobs()-rows(acor_grid),1,.))
st_store(.,"mean_LCI", mean_LCI\J(st_nobs()-rows(mean_LCI),1,.))
st_store(.,"acov_LCI", acov_LCI\J(st_nobs()-rows(acov_LCI),1,.))
st_store(.,"acor_LCI", acor_LCI\J(st_nobs()-rows(acor_LCI),1,.))
st_store(.,"mean_UCI", mean_UCI\J(st_nobs()-rows(mean_UCI),1,.))
st_store(.,"acov_UCI", acov_UCI\J(st_nobs()-rows(acov_UCI),1,.))
st_store(.,"acor_UCI", acor_UCI\J(st_nobs()-rows(acor_UCI),1,.))
}
// 2.3. Third-Order-Jackknife
/*computing TOJ empirical CDF estimate for T equivalent to 0 modulo 6
x point at which the empirical CDF is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ T/2
X22 vector of half-panel cross-sectional data based on time series (T/2 + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ T/3
X32 vector of one-third-panel cross-sectional data based on time series (T/3 + 1) ~ 2 * T/3
X33 vector of one-third-panel cross-sectional data based on time series 2 * T/3 + 1 ~ T
*/
function tojecdfest0(x, X, X21, X22, X31, X32, X33) {
// estimates
est = m_freq(x, X)
est21 = m_freq(x, X21)
est22 = m_freq(x, X22)
est31 = m_freq(x, X31)
est32 = m_freq(x, X32)
est33 = m_freq(x, X33)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33) / 3
// correction to ensure valid estimates
if (tojest < 0) {
tojest = 0
}
if (tojest > 1){
tojest = 1
}
return(tojest)
}
/* computing TOJ empirical CDF estimate for T equivalent to 1 modulo 6
x point at which the empricial CDF is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ floor(T/2)
X22 vector of half-panel cross-sectional data based on time series (floor(T/2) + 1) ~ T
X23 vector of half-panel cross-sectional data based on time series 1 ~ ceiling(T/2)
X24 vector of half-panel cross-sectional data based on time series (ceiling(T/2) + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X32 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3))
X33 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 1) ~ T
X34 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X35 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3) + 1)
X36 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 2) ~ T
X37 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X38 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * floor(T/3) + 1)
X39 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 2) ~ T
*/
function tojecdfest1(x, X, X21, X22, X23, X24, X31, X32, X33, X34, X35, X36, X37, X38, X39) {
// estimates
est = m_freq(x, X)
est21 = m_freq(x, X21)
est22 = m_freq(x, X22)
est23 = m_freq(x, X23)
est24 = m_freq(x, X24)
est31 = m_freq(x, X31)
est32 = m_freq(x, X32)
est33 = m_freq(x, X33)
est34 = m_freq(x, X34)
est35 = m_freq(x, X35)
est36 = m_freq(x, X36)
est37 = m_freq(x, X37)
est38 = m_freq(x, X38)
est39 = m_freq(x, X39)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
// correction to ensure valid estimates
if (tojest < 0) {
tojest = 0
}
if (tojest > 1){
tojest = 1
}
return(tojest)
}
/* computing TOJ empirical CDF estimate for T equivalent to 2 modulo 6
x point at which the empirical CDF is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ T/2
X22 vector of half-panel cross-sectional data based on time series (T/2 + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X32 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3) + 1)
X33 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3)) ~ T
X34 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X35 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * floor(T/3) + 1)
X36 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3)) ~ T
X37 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X38 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * ceiling(T/3))
X39 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3) + 1) ~ T
*/
function tojecdfest2(x, X, X21, X22, X31, X32, X33, X34, X35, X36, X37, X38, X39) {
// estimates
est = m_freq(x, X)
est21 = m_freq(x, X21)
est22 = m_freq(x, X22)
est31 = m_freq(x, X31)
est32 = m_freq(x, X32)
est33 = m_freq(x, X33)
est34 = m_freq(x, X34)
est35 = m_freq(x, X35)
est36 = m_freq(x, X36)
est37 = m_freq(x, X37)
est38 = m_freq(x, X38)
est39 = m_freq(x, X39)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
// correction to ensure valid estimates
if (tojest < 0) {
tojest = 0
}
if (tojest > 1){
tojest = 1
}
return(tojest)
}
/* computing TOJ empirical CDF estimate for T equivalent to 3 modulo 6
x point at which the empirical CDF is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ floor(T/2)
X22 vector of half-panel cross-sectional data based on time series (floor(T/2) + 1) ~ T
X23 vector of half-panel cross-sectional data based on time series 1 ~ ceiling(T/2)
X24 vector of half-panel cross-sectional data based on time series (ceiling(T/2) + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ T/3
X32 vector of one-third-panel cross-sectional data based on time series (T/3 + 1) ~ 2 * T/3
X33 vector of one-third-panel cross-sectional data based on time series 2 * T/3 + 1 ~ T
*/
function tojecdfest3(x, X, X21, X22, X23, X24, X31, X32, X33) {
// estimates
est = m_freq(x, X)
est21 = m_freq(x, X21)
est22 = m_freq(x, X22)
est23 = m_freq(x, X23)
est24 = m_freq(x, X24)
est31 = m_freq(x, X31)
est32 = m_freq(x, X32)
est33 = m_freq(x, X33)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33) / 3
// correction to ensure valid estimates
if (tojest < 0) {
tojest = 0
}
if (tojest > 1){
tojest = 1
}
return(tojest)
}
/* computing TOJ empirical CDF estimate for T equivalent to 4 modulo 6
x point at which the empirical CDF is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ T/2
X22 vector of half-panel cross-sectional data based on time series (T/2 + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X32 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3))
X33 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 1) ~ T
X34 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X35 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3) + 1)
X36 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 2) ~ T
X37 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X38 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * floor(T/3) + 1)
X39 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 2) ~ T
*/
function tojecdfest4(x, X, X21, X22, X31, X32, X33, X34, X35, X36, X37, X38, X39) {
// estimates
est = m_freq(x, X)
est21 = m_freq(x, X21)
est22 = m_freq(x, X22)
est31 = m_freq(x, X31)
est32 = m_freq(x, X32)
est33 = m_freq(x, X33)
est34 = m_freq(x, X34)
est35 = m_freq(x, X35)
est36 = m_freq(x, X36)
est37 = m_freq(x, X37)
est38 = m_freq(x, X38)
est39 = m_freq(x, X39)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
// correction to ensure valid estimates
if (tojest < 0) {
tojest = 0
}
if (tojest > 1){
tojest = 1
}
return(tojest)
}
/* computing TOJ empirical CDF estimate for T equivalent to 5 modulo 6
x point at which the empirical CDF is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ floor(T/2)
X22 vector of half-panel cross-sectional data based on time series (floor(T/2) + 1) ~ T
X23 vector of half-panel cross-sectional data based on time series 1 ~ ceiling(T/2)
X24 vector of half-panel cross-sectional data based on time series (ceiling(T/2) + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X32 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3) + 1)
X33 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3)) ~ T
X34 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X35 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * floor(T/3) + 1)
X36 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3)) ~ T
X37 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X38 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * ceiling(T/3))
X39 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3) + 1) ~ T
*/
function tojecdfest5(x, X, X21, X22, X23, X24, X31, X32, X33, X34, X35, X36, X37, X38, X39) {
// estimates
est = m_freq(x, X)
est21 = m_freq(x, X21)
est22 = m_freq(x, X22)
est23 = m_freq(x, X23)
est24 = m_freq(x, X24)
est31 = m_freq(x, X31)
est32 = m_freq(x, X32)
est33 = m_freq(x, X33)
est34 = m_freq(x, X34)
est35 = m_freq(x, X35)
est36 = m_freq(x, X36)
est37 = m_freq(x, X37)
est38 = m_freq(x, X38)
est39 = m_freq(x, X39)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
// correction to ensure valid estimates
if (tojest < 0) {
tojest = 0
}
if (tojest > 1){
tojest = 1
}
return(tojest)
}
function m_tojecdf(data, acov_order, acor_order, B) {
N = rows(data)
S = cols(data)
grid = 100
mean_est = J(N,1,0)
acov_est = J(N,1,0)
acor_est = J(N,1,0)
for (i=1; i<=N; i++){
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
}
if (mod(S,6) == 0) {
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::(S / 3)]
data32 = data[., (S / 3 + 1)::(2*S / 3)]
data33 = data[., (2 * S / 3 + 1)::S]
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i = 1; i <= grid; i++) {
mean_ecdf[i] = tojecdfest0(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est31, mean_est32, mean_est33)
acov_ecdf[i] = tojecdfest0(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est31, acov_est32, acov_est33)
acor_ecdf[i] = tojecdfest0(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est31, acor_est32, acor_est33)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = tojecdfest0(mean_grid[i], mean_est[index_boot], mean_est21[index_boot], mean_est22[index_boot], mean_est31[index_boot], mean_est32[index_boot], mean_est33[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = tojecdfest0(acov_grid[i], acov_est[index_boot], acov_est21[index_boot], acov_est22[index_boot], acov_est31[index_boot], acov_est32[index_boot], acov_est33[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = tojecdfest0(acor_grid[i], acor_est[index_boot], acor_est21[index_boot], acor_est22[index_boot], acor_est31[index_boot], acor_est32[index_boot], acor_est33[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
mean_data = sort((mean_ecdf, mean_LCI, mean_UCI), 1)
acov_data = sort((acov_ecdf, acov_LCI, acov_UCI), 1)
acor_data = sort((acor_ecdf, acor_LCI, acor_UCI), 1)
mean_ecdf = mean_data[,1]
mean_LCI = mean_data[,2]
mean_UCI = mean_data[,3]
acov_ecdf = acov_data[,1]
acov_LCI = acov_data[,2]
acov_UCI = acov_data[,3]
acor_ecdf = acor_data[,1]
acor_LCI = acor_data[,2]
acor_UCI = acor_data[,3]
} else if (mod(S,6)==1){
// split panel data for T equivalent to 1 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3))]
data33 = data[., (2 * floor(S / 3) + 1)::S]
data34 = data[., 1::floor(S / 3)]
data35 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * floor(S / 3) + 2)::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data39 = data[., (2 * floor(S / 3) + 2)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i = 1; i <= grid; i++) {
mean_ecdf[i] = tojecdfest1(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, mean_est36, mean_est37, mean_est38, mean_est39)
acov_ecdf[i] = tojecdfest1(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39)
acor_ecdf[i] = tojecdfest1(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, acor_est38, acor_est39)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = tojecdfest1(mean_grid[i], mean_est[index_boot], mean_est21[index_boot], mean_est22[index_boot], mean_est23[index_boot], mean_est24[index_boot], mean_est31[index_boot], mean_est32[index_boot], mean_est33[index_boot], mean_est34[index_boot], mean_est35[index_boot], mean_est36[index_boot], mean_est37[index_boot], mean_est38[index_boot], mean_est39[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = tojecdfest1(acov_grid[i], acov_est[index_boot], acov_est21[index_boot], acov_est22[index_boot], acov_est23[index_boot], acov_est24[index_boot], acov_est31[index_boot], acov_est32[index_boot], acov_est33[index_boot], acov_est34[index_boot], acov_est35[index_boot], acov_est36[index_boot], acov_est37[index_boot], acov_est38[index_boot], acov_est39[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = tojecdfest1(acor_grid[i], acor_est[index_boot], acor_est21[index_boot], acor_est22[index_boot], acor_est23[index_boot], acor_est24[index_boot], acor_est31[index_boot], acor_est32[index_boot], acor_est33[index_boot], acor_est34[index_boot], acor_est35[index_boot], acor_est36[index_boot], acor_est37[index_boot], acor_est38[index_boot], acor_est39[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
mean_data = sort((mean_ecdf, mean_LCI, mean_UCI), 1)
acov_data = sort((acov_ecdf, acov_LCI, acov_UCI), 1)
acor_data = sort((acor_ecdf, acor_LCI, acor_UCI), 1)
mean_ecdf = mean_data[,1]
mean_LCI = mean_data[,2]
mean_UCI = mean_data[,3]
acov_ecdf = acov_data[,1]
acov_LCI = acov_data[,2]
acov_UCI = acov_data[,3]
acor_ecdf = acor_data[,1]
acor_LCI = acor_data[,2]
acor_UCI = acor_data[,3]
} else if (mod(S,6)==2){
// split panel data for T equivalent to 2 modulo 6
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1) ]
data33 = data[., (2 * ceil(S / 3))::S]
data34 = data[., 1::ceil(S / 3)]
data35 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * ceil(S / 3))::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * ceil(S / 3))]
data39 = data[., (2 * ceil(S / 3) + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i = 1; i <= grid; i++) {
mean_ecdf[i] = tojecdfest2(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, mean_est36, mean_est37, mean_est38, mean_est39)
acov_ecdf[i] = tojecdfest2(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est31, acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39)
acor_ecdf[i] = tojecdfest2(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, acor_est38, acor_est39)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = tojecdfest2(mean_grid[i], mean_est[index_boot], mean_est21[index_boot], mean_est22[index_boot], mean_est31[index_boot], mean_est32[index_boot], mean_est33[index_boot], mean_est34[index_boot], mean_est35[index_boot], mean_est36[index_boot], mean_est37[index_boot], mean_est38[index_boot], mean_est39[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = tojecdfest2(acov_grid[i], acov_est[index_boot], acov_est21[index_boot], acov_est22[index_boot], acov_est31[index_boot], acov_est32[index_boot], acov_est33[index_boot], acov_est34[index_boot], acov_est35[index_boot], acov_est36[index_boot], acov_est37[index_boot], acov_est38[index_boot], acov_est39[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = tojecdfest2(acor_grid[i], acor_est[index_boot], acor_est21[index_boot], acor_est22[index_boot], acor_est31[index_boot], acor_est32[index_boot], acor_est33[index_boot], acor_est34[index_boot], acor_est35[index_boot], acor_est36[index_boot], acor_est37[index_boot], acor_est38[index_boot], acor_est39[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
mean_data = sort((mean_ecdf, mean_LCI, mean_UCI), 1)
acov_data = sort((acov_ecdf, acov_LCI, acov_UCI), 1)
acor_data = sort((acor_ecdf, acor_LCI, acor_UCI), 1)
mean_ecdf = mean_data[,1]
mean_LCI = mean_data[,2]
mean_UCI = mean_data[,3]
acov_ecdf = acov_data[,1]
acov_LCI = acov_data[,2]
acov_UCI = acov_data[,3]
acor_ecdf = acor_data[,1]
acor_LCI = acor_data[,2]
acor_UCI = acor_data[,3]
} else if (mod(S,6) == 3) {
// split panel data for T equivalent to 3 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::(S / 3)]
data32 = data[., (S / 3 + 1)::(2*S / 3)]
data33 = data[., (2 * S / 3 + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i = 1; i <= grid; i++) {
mean_ecdf[i] = tojecdfest3(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33)
acov_ecdf[i] = tojecdfest3(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, acov_est32, acov_est33)
acor_ecdf[i] = tojecdfest3(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = tojecdfest3(mean_grid[i], mean_est[index_boot], mean_est21[index_boot], mean_est22[index_boot], mean_est23[index_boot], mean_est24[index_boot], mean_est31[index_boot], mean_est32[index_boot], mean_est33[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = tojecdfest3(acov_grid[i], acov_est[index_boot], acov_est21[index_boot], acov_est22[index_boot], acov_est23[index_boot], acov_est24[index_boot], acov_est31[index_boot], acov_est32[index_boot], acov_est33[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = tojecdfest3(acor_grid[i], acor_est[index_boot], acor_est21[index_boot], acor_est22[index_boot], acor_est23[index_boot], acor_est24[index_boot], acor_est31[index_boot], acor_est32[index_boot], acor_est33[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
mean_data = sort((mean_ecdf, mean_LCI, mean_UCI), 1)
acov_data = sort((acov_ecdf, acov_LCI, acov_UCI), 1)
acor_data = sort((acor_ecdf, acor_LCI, acor_UCI), 1)
mean_ecdf = mean_data[,1]
mean_LCI = mean_data[,2]
mean_UCI = mean_data[,3]
acov_ecdf = acov_data[,1]
acov_LCI = acov_data[,2]
acov_UCI = acov_data[,3]
acor_ecdf = acor_data[,1]
acor_LCI = acor_data[,2]
acor_UCI = acor_data[,3]
} else if (mod(S,6)==4) {
// split panel data for T equivalent to 4 modulo 6
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3))]
data33 = data[., (2 * floor(S / 3) + 1)::S]
data34 = data[., 1::floor(S / 3)]
data35 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * floor(S / 3) + 2)::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data39 = data[., (2 * floor(S / 3) + 2)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i = 1; i <= grid; i++) {
mean_ecdf[i] = tojecdfest4(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, mean_est36, mean_est37, mean_est38, mean_est39)
acov_ecdf[i] = tojecdfest4(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est31, acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39)
acor_ecdf[i] = tojecdfest4(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, acor_est38, acor_est39)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = tojecdfest4(mean_grid[i], mean_est[index_boot], mean_est21[index_boot], mean_est22[index_boot], mean_est31[index_boot], mean_est32[index_boot], mean_est33[index_boot], mean_est34[index_boot], mean_est35[index_boot], mean_est36[index_boot], mean_est37[index_boot], mean_est38[index_boot], mean_est39[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = tojecdfest4(acov_grid[i], acov_est[index_boot], acov_est21[index_boot], acov_est22[index_boot], acov_est31[index_boot], acov_est32[index_boot], acov_est33[index_boot], acov_est34[index_boot], acov_est35[index_boot], acov_est36[index_boot], acov_est37[index_boot], acov_est38[index_boot], acov_est39[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = tojecdfest4(acor_grid[i], acor_est[index_boot], acor_est21[index_boot], acor_est22[index_boot], acor_est31[index_boot], acor_est32[index_boot], acor_est33[index_boot], acor_est34[index_boot], acor_est35[index_boot], acor_est36[index_boot], acor_est37[index_boot], acor_est38[index_boot], acor_est39[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
mean_data = sort((mean_ecdf, mean_LCI, mean_UCI), 1)
acov_data = sort((acov_ecdf, acov_LCI, acov_UCI), 1)
acor_data = sort((acor_ecdf, acor_LCI, acor_UCI), 1)
mean_ecdf = mean_data[,1]
mean_LCI = mean_data[,2]
mean_UCI = mean_data[,3]
acov_ecdf = acov_data[,1]
acov_LCI = acov_data[,2]
acov_UCI = acov_data[,3]
acor_ecdf = acor_data[,1]
acor_LCI = acor_data[,2]
acor_UCI = acor_data[,3]
} else {
// split panel data for T equivalent to 5 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1) ]
data33 = data[., (2 * ceil(S / 3))::S]
data34 = data[., 1::ceil(S / 3)]
data35 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * ceil(S / 3))::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * ceil(S / 3))]
data39 = data[., (2 * ceil(S / 3) + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_ecdf = J(grid, 1, 0)
acov_ecdf = J(grid, 1, 0)
acor_ecdf = J(grid, 1, 0)
mean_UCI = J(grid, 1, 0)
acov_UCI = J(grid, 1, 0)
acor_UCI = J(grid, 1, 0)
mean_LCI = J(grid, 1, 0)
acov_LCI = J(grid, 1, 0)
acor_LCI = J(grid, 1, 0)
for (i = 1; i <= grid; i++) {
mean_ecdf[i] = tojecdfest5(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, mean_est36, mean_est37, mean_est38, mean_est39)
acov_ecdf[i] = tojecdfest5(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39)
acor_ecdf[i] = tojecdfest5(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, acor_est38, acor_est39)
estimate_boot = J(B, 3, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,1] = tojecdfest5(mean_grid[i], mean_est[index_boot], mean_est21[index_boot], mean_est22[index_boot], mean_est23[index_boot], mean_est24[index_boot], mean_est31[index_boot], mean_est32[index_boot], mean_est33[index_boot], mean_est34[index_boot], mean_est35[index_boot], mean_est36[index_boot], mean_est37[index_boot], mean_est38[index_boot], mean_est39[index_boot]) - mean_ecdf[i]
estimate_boot[b,2] = tojecdfest5(acov_grid[i], acov_est[index_boot], acov_est21[index_boot], acov_est22[index_boot], acov_est23[index_boot], acov_est24[index_boot], acov_est31[index_boot], acov_est32[index_boot], acov_est33[index_boot], acov_est34[index_boot], acov_est35[index_boot], acov_est36[index_boot], acov_est37[index_boot], acov_est38[index_boot], acov_est39[index_boot]) - acov_ecdf[i]
estimate_boot[b,3] = tojecdfest5(acor_grid[i], acor_est[index_boot], acor_est21[index_boot], acor_est22[index_boot], acor_est23[index_boot], acor_est24[index_boot], acor_est31[index_boot], acor_est32[index_boot], acor_est33[index_boot], acor_est34[index_boot], acor_est35[index_boot], acor_est36[index_boot], acor_est37[index_boot], acor_est38[index_boot], acor_est39[index_boot]) - acor_ecdf[i]
}
mean_LCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.025)
mean_UCI[i] = mean_ecdf[i] + mm_quantile(estimate_boot[, 1], 1, 0.975)
acov_LCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.025)
acov_UCI[i] = acov_ecdf[i] + mm_quantile(estimate_boot[, 2], 1, 0.975)
acor_LCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.025)
acor_UCI[i] = acor_ecdf[i] + mm_quantile(estimate_boot[, 3], 1, 0.975)
}
mean_data = sort((mean_ecdf, mean_LCI, mean_UCI), 1)
acov_data = sort((acov_ecdf, acov_LCI, acov_UCI), 1)
acor_data = sort((acor_ecdf, acor_LCI, acor_UCI), 1)
mean_ecdf = mean_data[,1]
mean_LCI = mean_data[,2]
mean_UCI = mean_data[,3]
acov_ecdf = acov_data[,1]
acov_LCI = acov_data[,2]
acov_UCI = acov_data[,3]
acor_ecdf = acor_data[,1]
acor_LCI = acor_data[,2]
acor_UCI = acor_data[,3]
}
temp=st_addvar("double", "mean_ecdf")
temp=st_addvar("double", "acov_ecdf")
temp=st_addvar("double", "acor_ecdf")
temp=st_addvar("double", "mean_grid")
temp=st_addvar("double", "acov_grid")
temp=st_addvar("double", "acor_grid")
temp=st_addvar("double", "mean_LCI")
temp=st_addvar("double", "acov_LCI")
temp=st_addvar("double", "acor_LCI")
temp=st_addvar("double", "mean_UCI")
temp=st_addvar("double", "acov_UCI")
temp=st_addvar("double", "acor_UCI")
st_addobs(max((0,grid - st_nobs())))
st_store(.,"mean_ecdf", mean_ecdf\J(st_nobs()-rows(mean_ecdf),1,.))
st_store(.,"acov_ecdf", acov_ecdf\J(st_nobs()-rows(acov_ecdf),1,.))
st_store(.,"acor_ecdf", acor_ecdf\J(st_nobs()-rows(acor_ecdf),1,.))
st_store(.,"mean_grid", mean_grid\J(st_nobs()-rows(mean_grid),1,.))
st_store(.,"acov_grid", acov_grid\J(st_nobs()-rows(acov_grid),1,.))
st_store(.,"acor_grid", acor_grid\J(st_nobs()-rows(acor_grid),1,.))
st_store(.,"mean_LCI", mean_LCI\J(st_nobs()-rows(mean_LCI),1,.))
st_store(.,"acov_LCI", acov_LCI\J(st_nobs()-rows(acov_LCI),1,.))
st_store(.,"acor_LCI", acor_LCI\J(st_nobs()-rows(acor_LCI),1,.))
st_store(.,"mean_UCI", mean_UCI\J(st_nobs()-rows(mean_UCI),1,.))
st_store(.,"acov_UCI", acov_UCI\J(st_nobs()-rows(acov_UCI),1,.))
st_store(.,"acor_UCI", acor_UCI\J(st_nobs()-rows(acor_UCI),1,.))
}
/// 3. Moment Estimation
// 3.1. Naive Estimation
function momentest (quantity, indices) {
mean_est = quantity[indices, 1]
acov_est = quantity[indices, 2]
acor_est = quantity[indices, 3]
mean_mean = sum(mean_est) / length(mean_est)
acov_mean = sum(acov_est) / length(acov_est)
acor_mean = sum(acor_est) / length(acor_est)
mean_var = variance(mean_est)
acov_var = variance(acov_est)
acor_var = variance(acor_est)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
estimate = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
return(estimate)
}
function m_nemoment (data, acov_order, acor_order, B) {
N = rows(data)
S = cols(data)
level = 0.95
mean_est = J(N, 1, 0)
acov_est = J(N, 1, 0)
acor_est = J(N, 1, 0)
for (i = 1; i <= N; i++) {
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
}
equantity = (mean_est, acov_est, acor_est)
estimate_value = momentest(equantity, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,.] = momentest(equantity, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
st_matrix("est", estimate_value')
st_matrix("se", se)
st_matrix("ci", ci')
printf("\n")
printf("Estimates for Moments.\n")
printf("Parameters Estimate \n")
printf("________________________________________________________________ \n")
printf("Mean of Mean %f\n",estimate_value[1])
printf("Mean of Autocovariance %f\n",estimate_value[2])
printf("Mean of Autocorrelation %f\n",estimate_value[3])
printf("Variance of Mean %f\n",estimate_value[4])
printf("Variance of Autocovariance %f\n",estimate_value[5])
printf("Variance of Autocorrelation %f\n",estimate_value[6])
printf("Correlation between Mean and Autocovariance %f\n",estimate_value[7])
printf("Correlation between Mean and Autocorrelation %f\n",estimate_value[8])
printf("Correlation between Autocovariance and Autocorelation %f\n",estimate_value[9])
printf("\n")
printf("%f %% Confidence Intervals for Moments.\n", level*100)
printf("Parameters Low High\n")
printf("__________________________________________________________________________________\n")
printf("Mean of Mean %f %f\n",ci_1[1],ci_2[1])
printf("Mean of Autocovariance %f %f\n",ci_1[2],ci_2[2])
printf("Mean of Autocorrelation %f %f\n",ci_1[3],ci_2[3])
printf("Variance of Mean %f %f\n",ci_1[4],ci_2[4])
printf("Variance of Autocovariance %f %f\n",ci_1[5],ci_2[5])
printf("Variance of Autocorrelation %f %f\n",ci_1[6],ci_2[6])
printf("Correlation between Mean and Autocovariance %f %f\n",ci_1[7],ci_2[7])
printf("Correlation between Mean and Autocorrelation %f %f\n",ci_1[8],ci_2[8])
printf("Correlation between Autocovariance and Autocorelation %f %f\n",ci_1[9],ci_2[9])
printf("\n")
printf("Standard Errors for Moments.\n")
printf("Parameters Stanadard Errors \n")
printf("________________________________________________________________________\n")
printf("Mean of Mean %f\n",se[1])
printf("Mean of Autocovariance %f\n",se[2])
printf("Mean of Autocorrelation %f\n",se[3])
printf("Variance of Mean %f\n",se[4])
printf("Variance of Autocovariance %f\n",se[5])
printf("Variance of Autocorrelation %f\n",se[6])
printf("Correlation between Mean and Autocovariance %f\n",se[7])
printf("Correlation between Mean and Autocorrelation %f\n",se[8])
printf("Correlation between Autocovariance and Autocorelation %f\n",se[9])
}
// 3.2. Half-Panel-Jackknife Moment Estimation
function hpjmomentest1(quantity, indices){
mean_est = quantity[indices, 1]
mean_est1 = quantity[indices, 2]
mean_est2 = quantity[indices, 3]
acov_est = quantity[indices, 4]
acov_est1 = quantity[indices, 5]
acov_est2 = quantity[indices, 6]
acor_est = quantity[indices, 7]
acor_est1 = quantity[indices, 8]
acor_est2 = quantity[indices, 9]
mean_mean = sum(mean_est) / length(mean_est)
mean_mean1 = sum(mean_est1) / length(mean_est1)
mean_mean2 = sum(mean_est2) / length(mean_est2)
acov_mean = sum(acov_est) / length(acov_est)
acov_mean1 = sum(acov_est1) / length(acov_est1)
acov_mean2 = sum(acov_est2) / length(acov_est2)
acor_mean = sum(acor_est) / length(acor_est)
acor_mean1 = sum(acor_est1) / length(acor_est1)
acor_mean2 = sum(acor_est2) / length(acor_est2)
mean_var = variance(mean_est)
mean_var1 = variance(mean_est1)
mean_var2 = variance(mean_est2)
acov_var = variance(acov_est)
acov_var1 = variance(acov_est1)
acov_var2 = variance(acov_est2)
acor_var = variance(acor_est)
acor_var1 = variance(acor_est1)
acor_var2 = variance(acor_est2)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acov_cor1 = correlation((mean_est1, acov_est1))[2,1]
mean_acov_cor2 = correlation((mean_est2, acov_est2))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
mean_acor_cor1 = correlation((mean_est1, acor_est1))[2,1]
mean_acor_cor2 = correlation((mean_est2, acor_est2))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
acov_acor_cor1 = correlation((acov_est1, acor_est1))[2,1]
acov_acor_cor2 = correlation((acov_est2, acor_est2))[2,1]
estimate_1 = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
estimate_2 = (mean_mean1, acov_mean1, acor_mean1, mean_var1, acov_var1, acor_var1, mean_acov_cor1, mean_acor_cor1, acov_acor_cor1)
estimate_3 = (mean_mean2, acov_mean2, acor_mean2, mean_var2, acov_var2, acor_var2, mean_acov_cor2, mean_acor_cor2, acov_acor_cor2)
hpjestimate = 2 * estimate_1 - (estimate_2 + estimate_3) / 2
return(hpjestimate)
}
function hpjmomentest2(quantity, indices){
mean_est = quantity[indices, 1]
mean_est1 = quantity[indices, 2]
mean_est2 = quantity[indices, 3]
mean_est3 = quantity[indices, 4]
mean_est4 = quantity[indices, 5]
acov_est = quantity[indices, 6]
acov_est1 = quantity[indices, 7]
acov_est2 = quantity[indices, 8]
acov_est3 = quantity[indices, 9]
acov_est4 = quantity[indices, 10]
acor_est = quantity[indices, 11]
acor_est1 = quantity[indices, 12]
acor_est2 = quantity[indices, 13]
acor_est3 = quantity[indices, 14]
acor_est4 = quantity[indices, 15]
mean_mean = sum(mean_est)/length(mean_est)
mean_mean1 = sum(mean_est1)/length(mean_est1)
mean_mean2 = sum(mean_est2)/length(mean_est2)
mean_mean3 = sum(mean_est3)/length(mean_est3)
mean_mean4 = sum(mean_est4)/length(mean_est4)
acov_mean = sum(acov_est)/length(acov_est)
acov_mean1 = sum(acov_est1)/length(acov_est1)
acov_mean2 = sum(acov_est2)/length(acov_est2)
acov_mean3 = sum(acov_est3)/length(acov_est3)
acov_mean4 = sum(acov_est4)/length(acov_est4)
acor_mean = sum(acor_est)/length(acor_est)
acor_mean1 = sum(acor_est1)/length(acor_est1)
acor_mean2 = sum(acor_est2)/length(acor_est2)
acor_mean3 = sum(acor_est3)/length(acor_est3)
acor_mean4 = sum(acor_est4)/length(acor_est4)
mean_var = variance(mean_est)
mean_var1 = variance(mean_est1)
mean_var2 = variance(mean_est2)
mean_var3 = variance(mean_est3)
mean_var4 = variance(mean_est4)
acov_var = variance(acov_est)
acov_var1 = variance(acov_est1)
acov_var2 = variance(acov_est2)
acov_var3 = variance(acov_est3)
acov_var4 = variance(acov_est4)
acor_var = variance(acor_est)
acor_var1 = variance(acor_est1)
acor_var2 = variance(acor_est2)
acor_var3 = variance(acor_est3)
acor_var4 = variance(acor_est4)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acov_cor1 = correlation((mean_est1, acov_est1))[2,1]
mean_acov_cor2 = correlation((mean_est2, acov_est2))[2,1]
mean_acov_cor3 = correlation((mean_est3, acov_est3))[2,1]
mean_acov_cor4 = correlation((mean_est4, acov_est4))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
mean_acor_cor1 = correlation((mean_est1, acor_est1))[2,1]
mean_acor_cor2 = correlation((mean_est2, acor_est2))[2,1]
mean_acor_cor3 = correlation((mean_est3, acor_est3))[2,1]
mean_acor_cor4 = correlation((mean_est4, acor_est4))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
acov_acor_cor1 = correlation((acov_est1, acor_est1))[2,1]
acov_acor_cor2 = correlation((acov_est2, acor_est2))[2,1]
acov_acor_cor3 = correlation((acov_est3, acor_est3))[2,1]
acov_acor_cor4 = correlation((acov_est4, acor_est4))[2,1]
estimate_1 = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
estimate_2 = (mean_mean1, acov_mean1, acor_mean1, mean_var1, acov_var1, acor_var1, mean_acov_cor1, mean_acor_cor1, acov_acor_cor1)
estimate_3 = (mean_mean2, acov_mean2, acor_mean2, mean_var2, acov_var2, acor_var2, mean_acov_cor2, mean_acor_cor2, acov_acor_cor2)
estimate_4 = (mean_mean3, acov_mean3, acor_mean3, mean_var3, acov_var3, acor_var3, mean_acov_cor3, mean_acor_cor3, acov_acor_cor3)
estimate_5 = (mean_mean4, acov_mean4, acor_mean4, mean_var4, acov_var4, acor_var4, mean_acov_cor4, mean_acor_cor4, acov_acor_cor4)
hpjestimate = 2 * estimate_1 - (estimate_2 + estimate_3 + estimate_4 + estimate_5) / 4
return(hpjestimate)
}
function m_hpjmoment(data, acov_order, acor_order, B){
N = rows(data)
S = cols(data)
level = 0.95
if (mod(S,2)==0) {
mean_est = J(N, 1, 0)
acov_est = J(N, 1, 0)
acor_est = J(N, 1, 0)
data1 = data[,1::(S/2)]
data2 = data[,(S/2+1)::S]
mean_est1 = J(N, 1, 0)
mean_est2 = J(N, 1, 0)
acov_est1 = J(N, 1, 0)
acor_est1 = J(N, 1, 0)
acov_est2 = J(N, 1, 0)
acor_est2 = J(N, 1, 0)
for(i=1; i<=N ;i++){
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
mean_est1[i] = mean(data1[i,.]')
acov_est1[i] = mataacov(data1[i,.], acov_order)
acor_est1[i] = mataacor(data1[i,.], acor_order)
mean_est2[i] = mean(data2[i,.]')
acov_est2[i] = mataacov(data2[i,.], acov_order)
acor_est2[i] = mataacor(data2[i,.], acor_order)
}
equantity2 = (mean_est, mean_est1, mean_est2, acov_est, acov_est1, acov_est2, acor_est, acor_est1, acor_est2)
estimate_value = hpjmomentest1(equantity2, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,] = hpjmomentest1(equantity2, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
}
else{
mean_est = J(N, 1, 0)
acov_est = J(N, 1, 0)
acor_est = J(N, 1, 0)
data1 = data[., 1::floor(S/2)]
data2 = data[., (floor(S/2)+1)::S]
data3 = data[.,1::(floor(S/2)+1)]
data4 = data[.,(floor(S/2)+2)::S]
mean_est1 = J(N, 1, 0)
mean_est2 = J(N, 1, 0)
mean_est3 = J(N, 1, 0)
mean_est4 = J(N, 1, 0)
acov_est1 = J(N, 1, 0)
acov_est2 = J(N, 1, 0)
acov_est3 = J(N, 1, 0)
acov_est4 = J(N, 1, 0)
acor_est1 = J(N, 1, 0)
acor_est2 = J(N, 1, 0)
acor_est3 = J(N, 1, 0)
acor_est4 = J(N, 1, 0)
for(i=1; i<=N; i++){
mean_est[i] = mean(data[i,.]')
mean_est1[i] = mean(data1[i,.]')
mean_est2[i] = mean(data2[i,.]')
mean_est3[i] = mean(data3[i,.]')
mean_est4[i] = mean(data4[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acov_est1[i] = mataacov(data1[i,.], acov_order)
acov_est2[i] = mataacov(data2[i,.], acov_order)
acov_est3[i] = mataacov(data3[i,.], acov_order)
acov_est4[i] = mataacov(data4[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
acor_est1[i] = mataacor(data1[i,.], acor_order)
acor_est2[i] = mataacor(data2[i,.], acor_order)
acor_est3[i] = mataacor(data3[i,.], acor_order)
acor_est4[i] = mataacor(data4[i,.], acor_order)
}
equantity3 = (mean_est,mean_est1,mean_est2,mean_est3,mean_est4, ///
acov_est,acov_est1,acov_est2,acov_est3,acov_est4, ///
acor_est,acor_est1,acor_est2,acor_est3,acor_est4)
estimate_value = hpjmomentest2(equantity3, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,.] = hpjmomentest2(equantity3, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
}
st_matrix("est", estimate_value')
st_matrix("se", se)
st_matrix("ci", ci')
printf("\n")
printf("Estimates for Moments.\n")
printf("Parameters Estimate \n")
printf("________________________________________________________________ \n")
printf("Mean of Mean %f\n",estimate_value[1])
printf("Mean of Autocovariance %f\n",estimate_value[2])
printf("Mean of Autocorrelation %f\n",estimate_value[3])
printf("Variance of Mean %f\n",estimate_value[4])
printf("Variance of Autocovariance %f\n",estimate_value[5])
printf("Variance of Autocorrelation %f\n",estimate_value[6])
printf("Correlation between Mean and Autocovariance %f\n",estimate_value[7])
printf("Correlation between Mean and Autocorrelation %f\n",estimate_value[8])
printf("Correlation between Autocovariance and Autocorelation %f\n",estimate_value[9])
printf("\n")
printf("%f %% Confidence Intervals for Moments.\n", level*100)
printf("Parameters Low High\n")
printf("__________________________________________________________________________________\n")
printf("Mean of Mean %f %f\n",ci_1[1],ci_2[1])
printf("Mean of Autocovariance %f %f\n",ci_1[2],ci_2[2])
printf("Mean of Autocorrelation %f %f\n",ci_1[3],ci_2[3])
printf("Variance of Mean %f %f\n",ci_1[4],ci_2[4])
printf("Variance of Autocovariance %f %f\n",ci_1[5],ci_2[5])
printf("Variance of Autocorrelation %f %f\n",ci_1[6],ci_2[6])
printf("Correlation between Mean and Autocovariance %f %f\n",ci_1[7],ci_2[7])
printf("Correlation between Mean and Autocorrelation %f %f\n",ci_1[8],ci_2[8])
printf("Correlation between Autocovariance and Autocorelation %f %f\n",ci_1[9],ci_2[9])
printf("\n")
printf("Standard Errors for Moments.\n")
printf("Parameters Stanadard Errors \n")
printf("________________________________________________________________________\n")
printf("Mean of Mean %f\n",se[1])
printf("Mean of Autocovariance %f\n",se[2])
printf("Mean of Autocorrelation %f\n",se[3])
printf("Variance of Mean %f\n",se[4])
printf("Variance of Autocovariance %f\n",se[5])
printf("Variance of Autocorrelation %f\n",se[6])
printf("Correlation between Mean and Autocovariance %f\n",se[7])
printf("Correlation between Mean and Autocorrelation %f\n",se[8])
printf("Correlation between Autocovariance and Autocorelation %f\n",se[9])
}
// 3.3. Third-Order-Jackknife Moment Estimation
function tojmomentest0(quantity, indices){
mean_est = quantity[indices, 1]
mean_est21 = quantity[indices, 2]
mean_est22 = quantity[indices, 3]
mean_est31 = quantity[indices, 4]
mean_est32 = quantity[indices, 5]
mean_est33 = quantity[indices, 6]
acov_est = quantity[indices, 7]
acov_est21 = quantity[indices, 8]
acov_est22 = quantity[indices, 9]
acov_est31 = quantity[indices, 10]
acov_est32 = quantity[indices, 11]
acov_est33 = quantity[indices, 12]
acor_est = quantity[indices, 13]
acor_est21 = quantity[indices, 14]
acor_est22 = quantity[indices, 15]
acor_est31 = quantity[indices, 16]
acor_est32 = quantity[indices, 17]
acor_est33 = quantity[indices, 18]
mean_mean = sum(mean_est) / length(mean_est)
mean_mean21 = sum(mean_est21) / length(mean_est21)
mean_mean22 = sum(mean_est22) / length(mean_est22)
mean_mean31 = sum(mean_est31) / length(mean_est31)
mean_mean32 = sum(mean_est32) / length(mean_est32)
mean_mean33 = sum(mean_est33) / length(mean_est33)
acov_mean = sum(acov_est) / length(acov_est)
acov_mean21 = sum(acov_est21) / length(acov_est21)
acov_mean22 = sum(acov_est22) / length(acov_est22)
acov_mean31 = sum(acov_est31) / length(acov_est31)
acov_mean32 = sum(acov_est32) / length(acov_est32)
acov_mean33 = sum(acov_est33) / length(acov_est33)
acor_mean = sum(acor_est) / length(acor_est)
acor_mean21 = sum(acor_est21) / length(acor_est21)
acor_mean22 = sum(acor_est22) / length(acor_est22)
acor_mean31 = sum(acor_est31) / length(acor_est31)
acor_mean32 = sum(acor_est32) / length(acor_est32)
acor_mean33 = sum(acor_est33) / length(acor_est33)
mean_var = variance(mean_est)
mean_var21 = variance(mean_est21)
mean_var22 = variance(mean_est22)
mean_var31 = variance(mean_est31)
mean_var32 = variance(mean_est32)
mean_var33 = variance(mean_est33)
acov_var = variance(acov_est)
acov_var21 = variance(acov_est21)
acov_var22 = variance(acov_est22)
acov_var31 = variance(acov_est31)
acov_var32 = variance(acov_est32)
acov_var33 = variance(acov_est33)
acor_var = variance(acor_est)
acor_var21 = variance(acor_est21)
acor_var22 = variance(acor_est22)
acor_var31 = variance(acor_est31)
acor_var32 = variance(acor_est32)
acor_var33 = variance(acor_est33)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acov_cor21 = correlation((mean_est21, acov_est21))[2,1]
mean_acov_cor22 = correlation((mean_est22, acov_est22))[2,1]
mean_acov_cor31 = correlation((mean_est31, acov_est31))[2,1]
mean_acov_cor32 = correlation((mean_est32, acov_est32))[2,1]
mean_acov_cor33 = correlation((mean_est33, acov_est33))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
mean_acor_cor21 = correlation((mean_est21, acor_est21))[2,1]
mean_acor_cor22 = correlation((mean_est22, acor_est22))[2,1]
mean_acor_cor31 = correlation((mean_est31, acor_est31))[2,1]
mean_acor_cor32 = correlation((mean_est32, acor_est32))[2,1]
mean_acor_cor33 = correlation((mean_est33, acor_est33))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
acov_acor_cor21 = correlation((acov_est21, acor_est21))[2,1]
acov_acor_cor22 = correlation((acov_est22, acor_est22))[2,1]
acov_acor_cor31 = correlation((acov_est31, acor_est31))[2,1]
acov_acor_cor32 = correlation((acov_est32, acor_est32))[2,1]
acov_acor_cor33 = correlation((acov_est33, acor_est33))[2,1]
est = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
est21 = (mean_mean21, acov_mean21, acor_mean21, mean_var21, acov_var21, acor_var21, mean_acov_cor21, mean_acor_cor21, acov_acor_cor21)
est22 = (mean_mean22, acov_mean22, acor_mean22, mean_var22, acov_var22, acor_var22, mean_acov_cor22, mean_acor_cor22, acov_acor_cor22)
est31 = (mean_mean31, acov_mean31, acor_mean31, mean_var31, acov_var31, acor_var31, mean_acov_cor31, mean_acor_cor31, acov_acor_cor31)
est32 = (mean_mean32, acov_mean32, acor_mean32, mean_var32, acov_var32, acor_var32, mean_acov_cor32, mean_acor_cor32, acov_acor_cor32)
est33 = (mean_mean33, acov_mean33, acor_mean33, mean_var33, acov_var33, acor_var33, mean_acov_cor33, mean_acor_cor33, acov_acor_cor33)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33) / 3
return(tojest)
}
function tojmomentest1(quantity, indices){
mean_est = quantity[indices, 1]
mean_est21 = quantity[indices, 2]
mean_est22 = quantity[indices, 3]
mean_est23 = quantity[indices, 4]
mean_est24 = quantity[indices, 5]
mean_est31 = quantity[indices, 6]
mean_est32 = quantity[indices, 7]
mean_est33 = quantity[indices, 8]
mean_est34 = quantity[indices, 9]
mean_est35 = quantity[indices, 10]
mean_est36 = quantity[indices, 11]
mean_est37 = quantity[indices, 12]
mean_est38 = quantity[indices, 13]
mean_est39 = quantity[indices, 14]
acov_est = quantity[indices, 15]
acov_est21 = quantity[indices, 16]
acov_est22 = quantity[indices, 17]
acov_est23 = quantity[indices, 18]
acov_est24 = quantity[indices, 19]
acov_est31 = quantity[indices, 20]
acov_est32 = quantity[indices, 21]
acov_est33 = quantity[indices, 22]
acov_est34 = quantity[indices, 23]
acov_est35 = quantity[indices, 24]
acov_est36 = quantity[indices, 25]
acov_est37 = quantity[indices, 26]
acov_est38 = quantity[indices, 27]
acov_est39 = quantity[indices, 28]
acor_est = quantity[indices, 29]
acor_est21 = quantity[indices, 30]
acor_est22 = quantity[indices, 31]
acor_est23 = quantity[indices, 32]
acor_est24 = quantity[indices, 33]
acor_est31 = quantity[indices, 34]
acor_est32 = quantity[indices, 35]
acor_est33 = quantity[indices, 36]
acor_est34 = quantity[indices, 37]
acor_est35 = quantity[indices, 38]
acor_est36 = quantity[indices, 39]
acor_est37 = quantity[indices, 40]
acor_est38 = quantity[indices, 41]
acor_est39 = quantity[indices, 42]
mean_mean = sum(mean_est) / length(mean_est)
mean_mean21 = sum(mean_est21) / length(mean_est21)
mean_mean22 = sum(mean_est22) / length(mean_est22)
mean_mean23 = sum(mean_est23) / length(mean_est23)
mean_mean24 = sum(mean_est24) / length(mean_est24)
mean_mean31 = sum(mean_est31) / length(mean_est31)
mean_mean32 = sum(mean_est32) / length(mean_est32)
mean_mean33 = sum(mean_est33) / length(mean_est33)
mean_mean34 = sum(mean_est34) / length(mean_est34)
mean_mean35 = sum(mean_est35) / length(mean_est35)
mean_mean36 = sum(mean_est36) / length(mean_est36)
mean_mean37 = sum(mean_est37) / length(mean_est37)
mean_mean38 = sum(mean_est38) / length(mean_est38)
mean_mean39 = sum(mean_est39) / length(mean_est39)
acov_mean = sum(acov_est) / length(acov_est)
acov_mean21 = sum(acov_est21) / length(acov_est21)
acov_mean22 = sum(acov_est22) / length(acov_est22)
acov_mean23 = sum(acov_est23) / length(acov_est23)
acov_mean24 = sum(acov_est24) / length(acov_est24)
acov_mean31 = sum(acov_est31) / length(acov_est31)
acov_mean32 = sum(acov_est32) / length(acov_est32)
acov_mean33 = sum(acov_est33) / length(acov_est33)
acov_mean34 = sum(acov_est34) / length(acov_est34)
acov_mean35 = sum(acov_est35) / length(acov_est35)
acov_mean36 = sum(acov_est36) / length(acov_est36)
acov_mean37 = sum(acov_est37) / length(acov_est37)
acov_mean38 = sum(acov_est38) / length(acov_est38)
acov_mean39 = sum(acov_est39) / length(acov_est39)
acor_mean = sum(acor_est) / length(acor_est)
acor_mean21 = sum(acor_est21) / length(acor_est21)
acor_mean22 = sum(acor_est22) / length(acor_est22)
acor_mean23 = sum(acor_est23) / length(acor_est23)
acor_mean24 = sum(acor_est24) / length(acor_est24)
acor_mean31 = sum(acor_est31) / length(acor_est31)
acor_mean32 = sum(acor_est32) / length(acor_est32)
acor_mean33 = sum(acor_est33) / length(acor_est33)
acor_mean34 = sum(acor_est34) / length(acor_est34)
acor_mean35 = sum(acor_est35) / length(acor_est35)
acor_mean36 = sum(acor_est36) / length(acor_est36)
acor_mean37 = sum(acor_est37) / length(acor_est37)
acor_mean38 = sum(acor_est38) / length(acor_est38)
acor_mean39 = sum(acor_est39) / length(acor_est39)
mean_var = variance(mean_est)
mean_var21 = variance(mean_est21)
mean_var22 = variance(mean_est22)
mean_var23 = variance(mean_est23)
mean_var24 = variance(mean_est24)
mean_var31 = variance(mean_est31)
mean_var32 = variance(mean_est32)
mean_var33 = variance(mean_est33)
mean_var34 = variance(mean_est34)
mean_var35 = variance(mean_est35)
mean_var36 = variance(mean_est36)
mean_var37 = variance(mean_est37)
mean_var38 = variance(mean_est38)
mean_var39 = variance(mean_est39)
acov_var = variance(acov_est)
acov_var21 = variance(acov_est21)
acov_var22 = variance(acov_est22)
acov_var23 = variance(acov_est23)
acov_var24 = variance(acov_est24)
acov_var31 = variance(acov_est31)
acov_var32 = variance(acov_est32)
acov_var33 = variance(acov_est33)
acov_var34 = variance(acov_est34)
acov_var35 = variance(acov_est35)
acov_var36 = variance(acov_est36)
acov_var37 = variance(acov_est37)
acov_var38 = variance(acov_est38)
acov_var39 = variance(acov_est39)
acor_var = variance(acor_est)
acor_var21 = variance(acor_est21)
acor_var22 = variance(acor_est22)
acor_var23 = variance(acor_est23)
acor_var24 = variance(acor_est24)
acor_var31 = variance(acor_est31)
acor_var32 = variance(acor_est32)
acor_var33 = variance(acor_est33)
acor_var34 = variance(acor_est34)
acor_var35 = variance(acor_est35)
acor_var36 = variance(acor_est36)
acor_var37 = variance(acor_est37)
acor_var38 = variance(acor_est38)
acor_var39 = variance(acor_est39)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acov_cor21 = correlation((mean_est21, acov_est21))[2,1]
mean_acov_cor22 = correlation((mean_est22, acov_est22))[2,1]
mean_acov_cor23 = correlation((mean_est23, acov_est23))[2,1]
mean_acov_cor24 = correlation((mean_est24, acov_est24))[2,1]
mean_acov_cor31 = correlation((mean_est31, acov_est31))[2,1]
mean_acov_cor32 = correlation((mean_est32, acov_est32))[2,1]
mean_acov_cor33 = correlation((mean_est33, acov_est33))[2,1]
mean_acov_cor34 = correlation((mean_est34, acov_est34))[2,1]
mean_acov_cor35 = correlation((mean_est35, acov_est35))[2,1]
mean_acov_cor36 = correlation((mean_est36, acov_est36))[2,1]
mean_acov_cor37 = correlation((mean_est37, acov_est37))[2,1]
mean_acov_cor38 = correlation((mean_est38, acov_est38))[2,1]
mean_acov_cor39 = correlation((mean_est39, acov_est39))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
mean_acor_cor21 = correlation((mean_est21, acor_est21))[2,1]
mean_acor_cor22 = correlation((mean_est22, acor_est22))[2,1]
mean_acor_cor23 = correlation((mean_est23, acor_est23))[2,1]
mean_acor_cor24 = correlation((mean_est24, acor_est24))[2,1]
mean_acor_cor31 = correlation((mean_est31, acor_est31))[2,1]
mean_acor_cor32 = correlation((mean_est32, acor_est32))[2,1]
mean_acor_cor33 = correlation((mean_est33, acor_est33))[2,1]
mean_acor_cor34 = correlation((mean_est34, acor_est34))[2,1]
mean_acor_cor35 = correlation((mean_est35, acor_est35))[2,1]
mean_acor_cor36 = correlation((mean_est36, acor_est36))[2,1]
mean_acor_cor37 = correlation((mean_est37, acor_est37))[2,1]
mean_acor_cor38 = correlation((mean_est38, acor_est38))[2,1]
mean_acor_cor39 = correlation((mean_est39, acor_est39))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
acov_acor_cor21 = correlation((acov_est21, acor_est21))[2,1]
acov_acor_cor22 = correlation((acov_est22, acor_est22))[2,1]
acov_acor_cor23 = correlation((acov_est23, acor_est23))[2,1]
acov_acor_cor24 = correlation((acov_est24, acor_est24))[2,1]
acov_acor_cor31 = correlation((acov_est31, acor_est31))[2,1]
acov_acor_cor32 = correlation((acov_est32, acor_est32))[2,1]
acov_acor_cor33 = correlation((acov_est33, acor_est33))[2,1]
acov_acor_cor34 = correlation((acov_est34, acor_est34))[2,1]
acov_acor_cor35 = correlation((acov_est35, acor_est35))[2,1]
acov_acor_cor36 = correlation((acov_est36, acor_est36))[2,1]
acov_acor_cor37 = correlation((acov_est37, acor_est37))[2,1]
acov_acor_cor38 = correlation((acov_est38, acor_est38))[2,1]
acov_acor_cor39 = correlation((acov_est39, acor_est39))[2,1]
est = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
est21 = (mean_mean21, acov_mean21, acor_mean21, mean_var21, acov_var21, acor_var21, mean_acov_cor21, mean_acor_cor21, acov_acor_cor21)
est22 = (mean_mean22, acov_mean22, acor_mean22, mean_var22, acov_var22, acor_var22, mean_acov_cor22, mean_acor_cor22, acov_acor_cor22)
est23 = (mean_mean23, acov_mean23, acor_mean23, mean_var23, acov_var23, acor_var23, mean_acov_cor23, mean_acor_cor23, acov_acor_cor23)
est24 = (mean_mean24, acov_mean24, acor_mean24, mean_var24, acov_var24, acor_var24, mean_acov_cor24, mean_acor_cor24, acov_acor_cor24)
est31 = (mean_mean31, acov_mean31, acor_mean31, mean_var31, acov_var31, acor_var31, mean_acov_cor31, mean_acor_cor31, acov_acor_cor31)
est32 = (mean_mean32, acov_mean32, acor_mean32, mean_var32, acov_var32, acor_var32, mean_acov_cor32, mean_acor_cor32, acov_acor_cor32)
est33 = (mean_mean33, acov_mean33, acor_mean33, mean_var33, acov_var33, acor_var33, mean_acov_cor33, mean_acor_cor33, acov_acor_cor33)
est34 = (mean_mean34, acov_mean34, acor_mean34, mean_var34, acov_var34, acor_var34, mean_acov_cor34, mean_acor_cor34, acov_acor_cor34)
est35 = (mean_mean35, acov_mean35, acor_mean35, mean_var35, acov_var35, acor_var35, mean_acov_cor35, mean_acor_cor35, acov_acor_cor35)
est36 = (mean_mean36, acov_mean36, acor_mean36, mean_var36, acov_var36, acor_var36, mean_acov_cor36, mean_acor_cor36, acov_acor_cor36)
est37 = (mean_mean37, acov_mean37, acor_mean37, mean_var37, acov_var37, acor_var37, mean_acov_cor37, mean_acor_cor37, acov_acor_cor37)
est38 = (mean_mean38, acov_mean38, acor_mean38, mean_var38, acov_var38, acor_var38, mean_acov_cor38, mean_acor_cor38, acov_acor_cor38)
est39 = (mean_mean39, acov_mean39, acor_mean39, mean_var39, acov_var39, acor_var39, mean_acov_cor39, mean_acor_cor39, acov_acor_cor39)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
return(tojest)
}
function tojmomentest2(quantity, indices){
mean_est = quantity[indices, 1]
mean_est21 = quantity[indices, 2]
mean_est22 = quantity[indices, 3]
mean_est31 = quantity[indices, 4]
mean_est32 = quantity[indices, 5]
mean_est33 = quantity[indices, 6]
mean_est34 = quantity[indices, 7]
mean_est35 = quantity[indices, 8]
mean_est36 = quantity[indices, 9]
mean_est37 = quantity[indices, 10]
mean_est38 = quantity[indices, 11]
mean_est39 = quantity[indices, 12]
acov_est = quantity[indices, 13]
acov_est21 = quantity[indices, 14]
acov_est22 = quantity[indices, 15]
acov_est31 = quantity[indices, 16]
acov_est32 = quantity[indices, 17]
acov_est33 = quantity[indices, 18]
acov_est34 = quantity[indices, 19]
acov_est35 = quantity[indices, 20]
acov_est36 = quantity[indices, 21]
acov_est37 = quantity[indices, 22]
acov_est38 = quantity[indices, 23]
acov_est39 = quantity[indices, 24]
acor_est = quantity[indices, 25]
acor_est21 = quantity[indices, 26]
acor_est22 = quantity[indices, 27]
acor_est31 = quantity[indices, 28]
acor_est32 = quantity[indices, 29]
acor_est33 = quantity[indices, 30]
acor_est34 = quantity[indices, 31]
acor_est35 = quantity[indices, 32]
acor_est36 = quantity[indices, 33]
acor_est37 = quantity[indices, 34]
acor_est38 = quantity[indices, 35]
acor_est39 = quantity[indices, 36]
mean_mean = sum(mean_est) / length(mean_est)
mean_mean21 = sum(mean_est21) / length(mean_est21)
mean_mean22 = sum(mean_est22) / length(mean_est22)
mean_mean31 = sum(mean_est31) / length(mean_est31)
mean_mean32 = sum(mean_est32) / length(mean_est32)
mean_mean33 = sum(mean_est33) / length(mean_est33)
mean_mean34 = sum(mean_est34) / length(mean_est34)
mean_mean35 = sum(mean_est35) / length(mean_est35)
mean_mean36 = sum(mean_est36) / length(mean_est36)
mean_mean37 = sum(mean_est37) / length(mean_est37)
mean_mean38 = sum(mean_est38) / length(mean_est38)
mean_mean39 = sum(mean_est39) / length(mean_est39)
acov_mean = sum(acov_est) / length(acov_est)
acov_mean21 = sum(acov_est21) / length(acov_est21)
acov_mean22 = sum(acov_est22) / length(acov_est22)
acov_mean31 = sum(acov_est31) / length(acov_est31)
acov_mean32 = sum(acov_est32) / length(acov_est32)
acov_mean33 = sum(acov_est33) / length(acov_est33)
acov_mean34 = sum(acov_est34) / length(acov_est34)
acov_mean35 = sum(acov_est35) / length(acov_est35)
acov_mean36 = sum(acov_est36) / length(acov_est36)
acov_mean37 = sum(acov_est37) / length(acov_est37)
acov_mean38 = sum(acov_est38) / length(acov_est38)
acov_mean39 = sum(acov_est39) / length(acov_est39)
acor_mean = sum(acor_est) / length(acor_est)
acor_mean21 = sum(acor_est21) / length(acor_est21)
acor_mean22 = sum(acor_est22) / length(acor_est22)
acor_mean31 = sum(acor_est31) / length(acor_est31)
acor_mean32 = sum(acor_est32) / length(acor_est32)
acor_mean33 = sum(acor_est33) / length(acor_est33)
acor_mean34 = sum(acor_est34) / length(acor_est34)
acor_mean35 = sum(acor_est35) / length(acor_est35)
acor_mean36 = sum(acor_est36) / length(acor_est36)
acor_mean37 = sum(acor_est37) / length(acor_est37)
acor_mean38 = sum(acor_est38) / length(acor_est38)
acor_mean39 = sum(acor_est39) / length(acor_est39)
mean_var = variance(mean_est)
mean_var21 = variance(mean_est21)
mean_var22 = variance(mean_est22)
mean_var31 = variance(mean_est31)
mean_var32 = variance(mean_est32)
mean_var33 = variance(mean_est33)
mean_var34 = variance(mean_est34)
mean_var35 = variance(mean_est35)
mean_var36 = variance(mean_est36)
mean_var37 = variance(mean_est37)
mean_var38 = variance(mean_est38)
mean_var39 = variance(mean_est39)
acov_var = variance(acov_est)
acov_var21 = variance(acov_est21)
acov_var22 = variance(acov_est22)
acov_var31 = variance(acov_est31)
acov_var32 = variance(acov_est32)
acov_var33 = variance(acov_est33)
acov_var34 = variance(acov_est34)
acov_var35 = variance(acov_est35)
acov_var36 = variance(acov_est36)
acov_var37 = variance(acov_est37)
acov_var38 = variance(acov_est38)
acov_var39 = variance(acov_est39)
acor_var = variance(acor_est)
acor_var21 = variance(acor_est21)
acor_var22 = variance(acor_est22)
acor_var31 = variance(acor_est31)
acor_var32 = variance(acor_est32)
acor_var33 = variance(acor_est33)
acor_var34 = variance(acor_est34)
acor_var35 = variance(acor_est35)
acor_var36 = variance(acor_est36)
acor_var37 = variance(acor_est37)
acor_var38 = variance(acor_est38)
acor_var39 = variance(acor_est39)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acov_cor21 = correlation((mean_est21, acov_est21))[2,1]
mean_acov_cor22 = correlation((mean_est22, acov_est22))[2,1]
mean_acov_cor31 = correlation((mean_est31, acov_est31))[2,1]
mean_acov_cor32 = correlation((mean_est32, acov_est32))[2,1]
mean_acov_cor33 = correlation((mean_est33, acov_est33))[2,1]
mean_acov_cor34 = correlation((mean_est34, acov_est34))[2,1]
mean_acov_cor35 = correlation((mean_est35, acov_est35))[2,1]
mean_acov_cor36 = correlation((mean_est36, acov_est36))[2,1]
mean_acov_cor37 = correlation((mean_est37, acov_est37))[2,1]
mean_acov_cor38 = correlation((mean_est38, acov_est38))[2,1]
mean_acov_cor39 = correlation((mean_est39, acov_est39))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
mean_acor_cor21 = correlation((mean_est21, acor_est21))[2,1]
mean_acor_cor22 = correlation((mean_est22, acor_est22))[2,1]
mean_acor_cor31 = correlation((mean_est31, acor_est31))[2,1]
mean_acor_cor32 = correlation((mean_est32, acor_est32))[2,1]
mean_acor_cor33 = correlation((mean_est33, acor_est33))[2,1]
mean_acor_cor34 = correlation((mean_est34, acor_est34))[2,1]
mean_acor_cor35 = correlation((mean_est35, acor_est35))[2,1]
mean_acor_cor36 = correlation((mean_est36, acor_est36))[2,1]
mean_acor_cor37 = correlation((mean_est37, acor_est37))[2,1]
mean_acor_cor38 = correlation((mean_est38, acor_est38))[2,1]
mean_acor_cor39 = correlation((mean_est39, acor_est39))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
acov_acor_cor21 = correlation((acov_est21, acor_est21))[2,1]
acov_acor_cor22 = correlation((acov_est22, acor_est22))[2,1]
acov_acor_cor31 = correlation((acov_est31, acor_est31))[2,1]
acov_acor_cor32 = correlation((acov_est32, acor_est32))[2,1]
acov_acor_cor33 = correlation((acov_est33, acor_est33))[2,1]
acov_acor_cor34 = correlation((acov_est34, acor_est34))[2,1]
acov_acor_cor35 = correlation((acov_est35, acor_est35))[2,1]
acov_acor_cor36 = correlation((acov_est36, acor_est36))[2,1]
acov_acor_cor37 = correlation((acov_est37, acor_est37))[2,1]
acov_acor_cor38 = correlation((acov_est38, acor_est38))[2,1]
acov_acor_cor39 = correlation((acov_est39, acor_est39))[2,1]
est = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
est21 = (mean_mean21, acov_mean21, acor_mean21, mean_var21, acov_var21, acor_var21, mean_acov_cor21, mean_acor_cor21, acov_acor_cor21)
est22 = (mean_mean22, acov_mean22, acor_mean22, mean_var22, acov_var22, acor_var22, mean_acov_cor22, mean_acor_cor22, acov_acor_cor22)
est31 = (mean_mean31, acov_mean31, acor_mean31, mean_var31, acov_var31, acor_var31, mean_acov_cor31, mean_acor_cor31, acov_acor_cor31)
est32 = (mean_mean32, acov_mean32, acor_mean32, mean_var32, acov_var32, acor_var32, mean_acov_cor32, mean_acor_cor32, acov_acor_cor32)
est33 = (mean_mean33, acov_mean33, acor_mean33, mean_var33, acov_var33, acor_var33, mean_acov_cor33, mean_acor_cor33, acov_acor_cor33)
est34 = (mean_mean34, acov_mean34, acor_mean34, mean_var34, acov_var34, acor_var34, mean_acov_cor34, mean_acor_cor34, acov_acor_cor34)
est35 = (mean_mean35, acov_mean35, acor_mean35, mean_var35, acov_var35, acor_var35, mean_acov_cor35, mean_acor_cor35, acov_acor_cor35)
est36 = (mean_mean36, acov_mean36, acor_mean36, mean_var36, acov_var36, acor_var36, mean_acov_cor36, mean_acor_cor36, acov_acor_cor36)
est37 = (mean_mean37, acov_mean37, acor_mean37, mean_var37, acov_var37, acor_var37, mean_acov_cor37, mean_acor_cor37, acov_acor_cor37)
est38 = (mean_mean38, acov_mean38, acor_mean38, mean_var38, acov_var38, acor_var38, mean_acov_cor38, mean_acor_cor38, acov_acor_cor38)
est39 = (mean_mean39, acov_mean39, acor_mean39, mean_var39, acov_var39, acor_var39, mean_acov_cor39, mean_acor_cor39, acov_acor_cor39)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
return(tojest)
}
function tojmomentest3(quantity, indices){
mean_est = quantity[indices, 1]
mean_est21 = quantity[indices, 2]
mean_est22 = quantity[indices, 3]
mean_est23 = quantity[indices, 4]
mean_est24 = quantity[indices, 5]
mean_est31 = quantity[indices, 6]
mean_est32 = quantity[indices, 7]
mean_est33 = quantity[indices, 8]
acov_est = quantity[indices, 9]
acov_est21 = quantity[indices, 10]
acov_est22 = quantity[indices, 11]
acov_est23 = quantity[indices, 12]
acov_est24 = quantity[indices, 13]
acov_est31 = quantity[indices, 14]
acov_est32 = quantity[indices, 15]
acov_est33 = quantity[indices, 16]
acor_est = quantity[indices, 17]
acor_est21 = quantity[indices, 18]
acor_est22 = quantity[indices, 19]
acor_est23 = quantity[indices, 20]
acor_est24 = quantity[indices, 21]
acor_est31 = quantity[indices, 22]
acor_est32 = quantity[indices, 23]
acor_est33 = quantity[indices, 24]
mean_mean = sum(mean_est) / length(mean_est)
mean_mean21 = sum(mean_est21) / length(mean_est21)
mean_mean22 = sum(mean_est22) / length(mean_est22)
mean_mean23 = sum(mean_est23) / length(mean_est23)
mean_mean24 = sum(mean_est24) / length(mean_est24)
mean_mean31 = sum(mean_est31) / length(mean_est31)
mean_mean32 = sum(mean_est32) / length(mean_est32)
mean_mean33 = sum(mean_est33) / length(mean_est33)
acov_mean = sum(acov_est) / length(acov_est)
acov_mean21 = sum(acov_est21) / length(acov_est21)
acov_mean22 = sum(acov_est22) / length(acov_est22)
acov_mean23 = sum(acov_est23) / length(acov_est23)
acov_mean24 = sum(acov_est24) / length(acov_est24)
acov_mean31 = sum(acov_est31) / length(acov_est31)
acov_mean32 = sum(acov_est32) / length(acov_est32)
acov_mean33 = sum(acov_est33) / length(acov_est33)
acor_mean = sum(acor_est) / length(acor_est)
acor_mean21 = sum(acor_est21) / length(acor_est21)
acor_mean22 = sum(acor_est22) / length(acor_est22)
acor_mean23 = sum(acor_est23) / length(acor_est23)
acor_mean24 = sum(acor_est24) / length(acor_est24)
acor_mean31 = sum(acor_est31) / length(acor_est31)
acor_mean32 = sum(acor_est32) / length(acor_est32)
acor_mean33 = sum(acor_est33) / length(acor_est33)
mean_var = variance(mean_est)
mean_var21 = variance(mean_est21)
mean_var22 = variance(mean_est22)
mean_var23 = variance(mean_est23)
mean_var24 = variance(mean_est24)
mean_var31 = variance(mean_est31)
mean_var32 = variance(mean_est32)
mean_var33 = variance(mean_est33)
acov_var = variance(acov_est)
acov_var21 = variance(acov_est21)
acov_var22 = variance(acov_est22)
acov_var23 = variance(acov_est23)
acov_var24 = variance(acov_est24)
acov_var31 = variance(acov_est31)
acov_var32 = variance(acov_est32)
acov_var33 = variance(acov_est33)
acor_var = variance(acor_est)
acor_var21 = variance(acor_est21)
acor_var22 = variance(acor_est22)
acor_var23 = variance(acor_est23)
acor_var24 = variance(acor_est24)
acor_var31 = variance(acor_est31)
acor_var32 = variance(acor_est32)
acor_var33 = variance(acor_est33)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acov_cor21 = correlation((mean_est21, acov_est21))[2,1]
mean_acov_cor22 = correlation((mean_est22, acov_est22))[2,1]
mean_acov_cor23 = correlation((mean_est23, acov_est23))[2,1]
mean_acov_cor24 = correlation((mean_est24, acov_est24))[2,1]
mean_acov_cor31 = correlation((mean_est31, acov_est31))[2,1]
mean_acov_cor32 = correlation((mean_est32, acov_est32))[2,1]
mean_acov_cor33 = correlation((mean_est33, acov_est33))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
mean_acor_cor21 = correlation((mean_est21, acor_est21))[2,1]
mean_acor_cor22 = correlation((mean_est22, acor_est22))[2,1]
mean_acor_cor23 = correlation((mean_est23, acor_est23))[2,1]
mean_acor_cor24 = correlation((mean_est24, acor_est24))[2,1]
mean_acor_cor31 = correlation((mean_est31, acor_est31))[2,1]
mean_acor_cor32 = correlation((mean_est32, acor_est32))[2,1]
mean_acor_cor33 = correlation((mean_est33, acor_est33))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
acov_acor_cor21 = correlation((acov_est21, acor_est21))[2,1]
acov_acor_cor22 = correlation((acov_est22, acor_est22))[2,1]
acov_acor_cor23 = correlation((acov_est23, acor_est23))[2,1]
acov_acor_cor24 = correlation((acov_est24, acor_est24))[2,1]
acov_acor_cor31 = correlation((acov_est31, acor_est31))[2,1]
acov_acor_cor32 = correlation((acov_est32, acor_est32))[2,1]
acov_acor_cor33 = correlation((acov_est33, acor_est33))[2,1]
est = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
est21 = (mean_mean21, acov_mean21, acor_mean21, mean_var21, acov_var21, acor_var21, mean_acov_cor21, mean_acor_cor21, acov_acor_cor21)
est22 = (mean_mean22, acov_mean22, acor_mean22, mean_var22, acov_var22, acor_var22, mean_acov_cor22, mean_acor_cor22, acov_acor_cor22)
est23 = (mean_mean23, acov_mean23, acor_mean23, mean_var23, acov_var23, acor_var23, mean_acov_cor23, mean_acor_cor23, acov_acor_cor23)
est24 = (mean_mean24, acov_mean24, acor_mean24, mean_var24, acov_var24, acor_var24, mean_acov_cor24, mean_acor_cor24, acov_acor_cor24)
est31 = (mean_mean31, acov_mean31, acor_mean31, mean_var31, acov_var31, acor_var31, mean_acov_cor31, mean_acor_cor31, acov_acor_cor31)
est32 = (mean_mean32, acov_mean32, acor_mean32, mean_var32, acov_var32, acor_var32, mean_acov_cor32, mean_acor_cor32, acov_acor_cor32)
est33 = (mean_mean33, acov_mean33, acor_mean33, mean_var33, acov_var33, acor_var33, mean_acov_cor33, mean_acor_cor33, acov_acor_cor33)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33) / 3
return(tojest)
}
function tojmomentest4(quantity, indices){
mean_est = quantity[indices, 1]
mean_est21 = quantity[indices, 2]
mean_est22 = quantity[indices, 3]
mean_est31 = quantity[indices, 4]
mean_est32 = quantity[indices, 5]
mean_est33 = quantity[indices, 6]
mean_est34 = quantity[indices, 7]
mean_est35 = quantity[indices, 8]
mean_est36 = quantity[indices, 9]
mean_est37 = quantity[indices, 10]
mean_est38 = quantity[indices, 11]
mean_est39 = quantity[indices, 12]
acov_est = quantity[indices, 13]
acov_est21 = quantity[indices, 14]
acov_est22 = quantity[indices, 15]
acov_est31 = quantity[indices, 16]
acov_est32 = quantity[indices, 17]
acov_est33 = quantity[indices, 18]
acov_est34 = quantity[indices, 19]
acov_est35 = quantity[indices, 20]
acov_est36 = quantity[indices, 21]
acov_est37 = quantity[indices, 22]
acov_est38 = quantity[indices, 23]
acov_est39 = quantity[indices, 24]
acor_est = quantity[indices, 25]
acor_est21 = quantity[indices, 26]
acor_est22 = quantity[indices, 27]
acor_est31 = quantity[indices, 28]
acor_est32 = quantity[indices, 29]
acor_est33 = quantity[indices, 30]
acor_est34 = quantity[indices, 31]
acor_est35 = quantity[indices, 32]
acor_est36 = quantity[indices, 33]
acor_est37 = quantity[indices, 34]
acor_est38 = quantity[indices, 35]
acor_est39 = quantity[indices, 36]
mean_mean = sum(mean_est) / length(mean_est)
mean_mean21 = sum(mean_est21) / length(mean_est21)
mean_mean22 = sum(mean_est22) / length(mean_est22)
mean_mean31 = sum(mean_est31) / length(mean_est31)
mean_mean32 = sum(mean_est32) / length(mean_est32)
mean_mean33 = sum(mean_est33) / length(mean_est33)
mean_mean34 = sum(mean_est34) / length(mean_est34)
mean_mean35 = sum(mean_est35) / length(mean_est35)
mean_mean36 = sum(mean_est36) / length(mean_est36)
mean_mean37 = sum(mean_est37) / length(mean_est37)
mean_mean38 = sum(mean_est38) / length(mean_est38)
mean_mean39 = sum(mean_est39) / length(mean_est39)
acov_mean = sum(acov_est) / length(acov_est)
acov_mean21 = sum(acov_est21) / length(acov_est21)
acov_mean22 = sum(acov_est22) / length(acov_est22)
acov_mean31 = sum(acov_est31) / length(acov_est31)
acov_mean32 = sum(acov_est32) / length(acov_est32)
acov_mean33 = sum(acov_est33) / length(acov_est33)
acov_mean34 = sum(acov_est34) / length(acov_est34)
acov_mean35 = sum(acov_est35) / length(acov_est35)
acov_mean36 = sum(acov_est36) / length(acov_est36)
acov_mean37 = sum(acov_est37) / length(acov_est37)
acov_mean38 = sum(acov_est38) / length(acov_est38)
acov_mean39 = sum(acov_est39) / length(acov_est39)
acor_mean = sum(acor_est) / length(acor_est)
acor_mean21 = sum(acor_est21) / length(acor_est21)
acor_mean22 = sum(acor_est22) / length(acor_est22)
acor_mean31 = sum(acor_est31) / length(acor_est31)
acor_mean32 = sum(acor_est32) / length(acor_est32)
acor_mean33 = sum(acor_est33) / length(acor_est33)
acor_mean34 = sum(acor_est34) / length(acor_est34)
acor_mean35 = sum(acor_est35) / length(acor_est35)
acor_mean36 = sum(acor_est36) / length(acor_est36)
acor_mean37 = sum(acor_est37) / length(acor_est37)
acor_mean38 = sum(acor_est38) / length(acor_est38)
acor_mean39 = sum(acor_est39) / length(acor_est39)
mean_var = variance(mean_est)
mean_var21 = variance(mean_est21)
mean_var22 = variance(mean_est22)
mean_var31 = variance(mean_est31)
mean_var32 = variance(mean_est32)
mean_var33 = variance(mean_est33)
mean_var34 = variance(mean_est34)
mean_var35 = variance(mean_est35)
mean_var36 = variance(mean_est36)
mean_var37 = variance(mean_est37)
mean_var38 = variance(mean_est38)
mean_var39 = variance(mean_est39)
acov_var = variance(acov_est)
acov_var21 = variance(acov_est21)
acov_var22 = variance(acov_est22)
acov_var31 = variance(acov_est31)
acov_var32 = variance(acov_est32)
acov_var33 = variance(acov_est33)
acov_var34 = variance(acov_est34)
acov_var35 = variance(acov_est35)
acov_var36 = variance(acov_est36)
acov_var37 = variance(acov_est37)
acov_var38 = variance(acov_est38)
acov_var39 = variance(acov_est39)
acor_var = variance(acor_est)
acor_var21 = variance(acor_est21)
acor_var22 = variance(acor_est22)
acor_var31 = variance(acor_est31)
acor_var32 = variance(acor_est32)
acor_var33 = variance(acor_est33)
acor_var34 = variance(acor_est34)
acor_var35 = variance(acor_est35)
acor_var36 = variance(acor_est36)
acor_var37 = variance(acor_est37)
acor_var38 = variance(acor_est38)
acor_var39 = variance(acor_est39)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acov_cor21 = correlation((mean_est21, acov_est21))[2,1]
mean_acov_cor22 = correlation((mean_est22, acov_est22))[2,1]
mean_acov_cor31 = correlation((mean_est31, acov_est31))[2,1]
mean_acov_cor32 = correlation((mean_est32, acov_est32))[2,1]
mean_acov_cor33 = correlation((mean_est33, acov_est33))[2,1]
mean_acov_cor34 = correlation((mean_est34, acov_est34))[2,1]
mean_acov_cor35 = correlation((mean_est35, acov_est35))[2,1]
mean_acov_cor36 = correlation((mean_est36, acov_est36))[2,1]
mean_acov_cor37 = correlation((mean_est37, acov_est37))[2,1]
mean_acov_cor38 = correlation((mean_est38, acov_est38))[2,1]
mean_acov_cor39 = correlation((mean_est39, acov_est39))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
mean_acor_cor21 = correlation((mean_est21, acor_est21))[2,1]
mean_acor_cor22 = correlation((mean_est22, acor_est22))[2,1]
mean_acor_cor31 = correlation((mean_est31, acor_est31))[2,1]
mean_acor_cor32 = correlation((mean_est32, acor_est32))[2,1]
mean_acor_cor33 = correlation((mean_est33, acor_est33))[2,1]
mean_acor_cor34 = correlation((mean_est34, acor_est34))[2,1]
mean_acor_cor35 = correlation((mean_est35, acor_est35))[2,1]
mean_acor_cor36 = correlation((mean_est36, acor_est36))[2,1]
mean_acor_cor37 = correlation((mean_est37, acor_est37))[2,1]
mean_acor_cor38 = correlation((mean_est38, acor_est38))[2,1]
mean_acor_cor39 = correlation((mean_est39, acor_est39))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
acov_acor_cor21 = correlation((acov_est21, acor_est21))[2,1]
acov_acor_cor22 = correlation((acov_est22, acor_est22))[2,1]
acov_acor_cor31 = correlation((acov_est31, acor_est31))[2,1]
acov_acor_cor32 = correlation((acov_est32, acor_est32))[2,1]
acov_acor_cor33 = correlation((acov_est33, acor_est33))[2,1]
acov_acor_cor34 = correlation((acov_est34, acor_est34))[2,1]
acov_acor_cor35 = correlation((acov_est35, acor_est35))[2,1]
acov_acor_cor36 = correlation((acov_est36, acor_est36))[2,1]
acov_acor_cor37 = correlation((acov_est37, acor_est37))[2,1]
acov_acor_cor38 = correlation((acov_est38, acor_est38))[2,1]
acov_acor_cor39 = correlation((acov_est39, acor_est39))[2,1]
est = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
est21 = (mean_mean21, acov_mean21, acor_mean21, mean_var21, acov_var21, acor_var21, mean_acov_cor21, mean_acor_cor21, acov_acor_cor21)
est22 = (mean_mean22, acov_mean22, acor_mean22, mean_var22, acov_var22, acor_var22, mean_acov_cor22, mean_acor_cor22, acov_acor_cor22)
est31 = (mean_mean31, acov_mean31, acor_mean31, mean_var31, acov_var31, acor_var31, mean_acov_cor31, mean_acor_cor31, acov_acor_cor31)
est32 = (mean_mean32, acov_mean32, acor_mean32, mean_var32, acov_var32, acor_var32, mean_acov_cor32, mean_acor_cor32, acov_acor_cor32)
est33 = (mean_mean33, acov_mean33, acor_mean33, mean_var33, acov_var33, acor_var33, mean_acov_cor33, mean_acor_cor33, acov_acor_cor33)
est34 = (mean_mean34, acov_mean34, acor_mean34, mean_var34, acov_var34, acor_var34, mean_acov_cor34, mean_acor_cor34, acov_acor_cor34)
est35 = (mean_mean35, acov_mean35, acor_mean35, mean_var35, acov_var35, acor_var35, mean_acov_cor35, mean_acor_cor35, acov_acor_cor35)
est36 = (mean_mean36, acov_mean36, acor_mean36, mean_var36, acov_var36, acor_var36, mean_acov_cor36, mean_acor_cor36, acov_acor_cor36)
est37 = (mean_mean37, acov_mean37, acor_mean37, mean_var37, acov_var37, acor_var37, mean_acov_cor37, mean_acor_cor37, acov_acor_cor37)
est38 = (mean_mean38, acov_mean38, acor_mean38, mean_var38, acov_var38, acor_var38, mean_acov_cor38, mean_acor_cor38, acov_acor_cor38)
est39 = (mean_mean39, acov_mean39, acor_mean39, mean_var39, acov_var39, acor_var39, mean_acov_cor39, mean_acor_cor39, acov_acor_cor39)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
return(tojest)
}
function tojmomentest5(quantity, indices){
mean_est = quantity[indices, 1]
mean_est21 = quantity[indices, 2]
mean_est22 = quantity[indices, 3]
mean_est23 = quantity[indices, 4]
mean_est24 = quantity[indices, 5]
mean_est31 = quantity[indices, 6]
mean_est32 = quantity[indices, 7]
mean_est33 = quantity[indices, 8]
mean_est34 = quantity[indices, 9]
mean_est35 = quantity[indices, 10]
mean_est36 = quantity[indices, 11]
mean_est37 = quantity[indices, 12]
mean_est38 = quantity[indices, 13]
mean_est39 = quantity[indices, 14]
acov_est = quantity[indices, 15]
acov_est21 = quantity[indices, 16]
acov_est22 = quantity[indices, 17]
acov_est23 = quantity[indices, 18]
acov_est24 = quantity[indices, 19]
acov_est31 = quantity[indices, 20]
acov_est32 = quantity[indices, 21]
acov_est33 = quantity[indices, 22]
acov_est34 = quantity[indices, 23]
acov_est35 = quantity[indices, 24]
acov_est36 = quantity[indices, 25]
acov_est37 = quantity[indices, 26]
acov_est38 = quantity[indices, 27]
acov_est39 = quantity[indices, 28]
acor_est = quantity[indices, 29]
acor_est21 = quantity[indices, 30]
acor_est22 = quantity[indices, 31]
acor_est23 = quantity[indices, 32]
acor_est24 = quantity[indices, 33]
acor_est31 = quantity[indices, 34]
acor_est32 = quantity[indices, 35]
acor_est33 = quantity[indices, 36]
acor_est34 = quantity[indices, 37]
acor_est35 = quantity[indices, 38]
acor_est36 = quantity[indices, 39]
acor_est37 = quantity[indices, 40]
acor_est38 = quantity[indices, 41]
acor_est39 = quantity[indices, 42]
mean_mean = sum(mean_est) / length(mean_est)
mean_mean21 = sum(mean_est21) / length(mean_est21)
mean_mean22 = sum(mean_est22) / length(mean_est22)
mean_mean23 = sum(mean_est23) / length(mean_est23)
mean_mean24 = sum(mean_est24) / length(mean_est24)
mean_mean31 = sum(mean_est31) / length(mean_est31)
mean_mean32 = sum(mean_est32) / length(mean_est32)
mean_mean33 = sum(mean_est33) / length(mean_est33)
mean_mean34 = sum(mean_est34) / length(mean_est34)
mean_mean35 = sum(mean_est35) / length(mean_est35)
mean_mean36 = sum(mean_est36) / length(mean_est36)
mean_mean37 = sum(mean_est37) / length(mean_est37)
mean_mean38 = sum(mean_est38) / length(mean_est38)
mean_mean39 = sum(mean_est39) / length(mean_est39)
acov_mean = sum(acov_est) / length(acov_est)
acov_mean21 = sum(acov_est21) / length(acov_est21)
acov_mean22 = sum(acov_est22) / length(acov_est22)
acov_mean23 = sum(acov_est23) / length(acov_est23)
acov_mean24 = sum(acov_est24) / length(acov_est24)
acov_mean31 = sum(acov_est31) / length(acov_est31)
acov_mean32 = sum(acov_est32) / length(acov_est32)
acov_mean33 = sum(acov_est33) / length(acov_est33)
acov_mean34 = sum(acov_est34) / length(acov_est34)
acov_mean35 = sum(acov_est35) / length(acov_est35)
acov_mean36 = sum(acov_est36) / length(acov_est36)
acov_mean37 = sum(acov_est37) / length(acov_est37)
acov_mean38 = sum(acov_est38) / length(acov_est38)
acov_mean39 = sum(acov_est39) / length(acov_est39)
acor_mean = sum(acor_est) / length(acor_est)
acor_mean21 = sum(acor_est21) / length(acor_est21)
acor_mean22 = sum(acor_est22) / length(acor_est22)
acor_mean23 = sum(acor_est23) / length(acor_est23)
acor_mean24 = sum(acor_est24) / length(acor_est24)
acor_mean31 = sum(acor_est31) / length(acor_est31)
acor_mean32 = sum(acor_est32) / length(acor_est32)
acor_mean33 = sum(acor_est33) / length(acor_est33)
acor_mean34 = sum(acor_est34) / length(acor_est34)
acor_mean35 = sum(acor_est35) / length(acor_est35)
acor_mean36 = sum(acor_est36) / length(acor_est36)
acor_mean37 = sum(acor_est37) / length(acor_est37)
acor_mean38 = sum(acor_est38) / length(acor_est38)
acor_mean39 = sum(acor_est39) / length(acor_est39)
mean_var = variance(mean_est)
mean_var21 = variance(mean_est21)
mean_var22 = variance(mean_est22)
mean_var23 = variance(mean_est23)
mean_var24 = variance(mean_est24)
mean_var31 = variance(mean_est31)
mean_var32 = variance(mean_est32)
mean_var33 = variance(mean_est33)
mean_var34 = variance(mean_est34)
mean_var35 = variance(mean_est35)
mean_var36 = variance(mean_est36)
mean_var37 = variance(mean_est37)
mean_var38 = variance(mean_est38)
mean_var39 = variance(mean_est39)
acov_var = variance(acov_est)
acov_var21 = variance(acov_est21)
acov_var22 = variance(acov_est22)
acov_var23 = variance(acov_est23)
acov_var24 = variance(acov_est24)
acov_var31 = variance(acov_est31)
acov_var32 = variance(acov_est32)
acov_var33 = variance(acov_est33)
acov_var34 = variance(acov_est34)
acov_var35 = variance(acov_est35)
acov_var36 = variance(acov_est36)
acov_var37 = variance(acov_est37)
acov_var38 = variance(acov_est38)
acov_var39 = variance(acov_est39)
acor_var = variance(acor_est)
acor_var21 = variance(acor_est21)
acor_var22 = variance(acor_est22)
acor_var23 = variance(acor_est23)
acor_var24 = variance(acor_est24)
acor_var31 = variance(acor_est31)
acor_var32 = variance(acor_est32)
acor_var33 = variance(acor_est33)
acor_var34 = variance(acor_est34)
acor_var35 = variance(acor_est35)
acor_var36 = variance(acor_est36)
acor_var37 = variance(acor_est37)
acor_var38 = variance(acor_est38)
acor_var39 = variance(acor_est39)
mean_acov_cor = correlation((mean_est, acov_est))[2,1]
mean_acov_cor21 = correlation((mean_est21, acov_est21))[2,1]
mean_acov_cor22 = correlation((mean_est22, acov_est22))[2,1]
mean_acov_cor23 = correlation((mean_est23, acov_est23))[2,1]
mean_acov_cor24 = correlation((mean_est24, acov_est24))[2,1]
mean_acov_cor31 = correlation((mean_est31, acov_est31))[2,1]
mean_acov_cor32 = correlation((mean_est32, acov_est32))[2,1]
mean_acov_cor33 = correlation((mean_est33, acov_est33))[2,1]
mean_acov_cor34 = correlation((mean_est34, acov_est34))[2,1]
mean_acov_cor35 = correlation((mean_est35, acov_est35))[2,1]
mean_acov_cor36 = correlation((mean_est36, acov_est36))[2,1]
mean_acov_cor37 = correlation((mean_est37, acov_est37))[2,1]
mean_acov_cor38 = correlation((mean_est38, acov_est38))[2,1]
mean_acov_cor39 = correlation((mean_est39, acov_est39))[2,1]
mean_acor_cor = correlation((mean_est, acor_est))[2,1]
mean_acor_cor21 = correlation((mean_est21, acor_est21))[2,1]
mean_acor_cor22 = correlation((mean_est22, acor_est22))[2,1]
mean_acor_cor23 = correlation((mean_est23, acor_est23))[2,1]
mean_acor_cor24 = correlation((mean_est24, acor_est24))[2,1]
mean_acor_cor31 = correlation((mean_est31, acor_est31))[2,1]
mean_acor_cor32 = correlation((mean_est32, acor_est32))[2,1]
mean_acor_cor33 = correlation((mean_est33, acor_est33))[2,1]
mean_acor_cor34 = correlation((mean_est34, acor_est34))[2,1]
mean_acor_cor35 = correlation((mean_est35, acor_est35))[2,1]
mean_acor_cor36 = correlation((mean_est36, acor_est36))[2,1]
mean_acor_cor37 = correlation((mean_est37, acor_est37))[2,1]
mean_acor_cor38 = correlation((mean_est38, acor_est38))[2,1]
mean_acor_cor39 = correlation((mean_est39, acor_est39))[2,1]
acov_acor_cor = correlation((acov_est, acor_est))[2,1]
acov_acor_cor21 = correlation((acov_est21, acor_est21))[2,1]
acov_acor_cor22 = correlation((acov_est22, acor_est22))[2,1]
acov_acor_cor23 = correlation((acov_est23, acor_est23))[2,1]
acov_acor_cor24 = correlation((acov_est24, acor_est24))[2,1]
acov_acor_cor31 = correlation((acov_est31, acor_est31))[2,1]
acov_acor_cor32 = correlation((acov_est32, acor_est32))[2,1]
acov_acor_cor33 = correlation((acov_est33, acor_est33))[2,1]
acov_acor_cor34 = correlation((acov_est34, acor_est34))[2,1]
acov_acor_cor35 = correlation((acov_est35, acor_est35))[2,1]
acov_acor_cor36 = correlation((acov_est36, acor_est36))[2,1]
acov_acor_cor37 = correlation((acov_est37, acor_est37))[2,1]
acov_acor_cor38 = correlation((acov_est38, acor_est38))[2,1]
acov_acor_cor39 = correlation((acov_est39, acor_est39))[2,1]
est = (mean_mean, acov_mean, acor_mean, mean_var, acov_var, acor_var, mean_acov_cor, mean_acor_cor, acov_acor_cor)
est21 = (mean_mean21, acov_mean21, acor_mean21, mean_var21, acov_var21, acor_var21, mean_acov_cor21, mean_acor_cor21, acov_acor_cor21)
est22 = (mean_mean22, acov_mean22, acor_mean22, mean_var22, acov_var22, acor_var22, mean_acov_cor22, mean_acor_cor22, acov_acor_cor22)
est23 = (mean_mean23, acov_mean23, acor_mean23, mean_var23, acov_var23, acor_var23, mean_acov_cor23, mean_acor_cor23, acov_acor_cor23)
est24 = (mean_mean24, acov_mean24, acor_mean24, mean_var24, acov_var24, acor_var24, mean_acov_cor24, mean_acor_cor24, acov_acor_cor24)
est31 = (mean_mean31, acov_mean31, acor_mean31, mean_var31, acov_var31, acor_var31, mean_acov_cor31, mean_acor_cor31, acov_acor_cor31)
est32 = (mean_mean32, acov_mean32, acor_mean32, mean_var32, acov_var32, acor_var32, mean_acov_cor32, mean_acor_cor32, acov_acor_cor32)
est33 = (mean_mean33, acov_mean33, acor_mean33, mean_var33, acov_var33, acor_var33, mean_acov_cor33, mean_acor_cor33, acov_acor_cor33)
est34 = (mean_mean34, acov_mean34, acor_mean34, mean_var34, acov_var34, acor_var34, mean_acov_cor34, mean_acor_cor34, acov_acor_cor34)
est35 = (mean_mean35, acov_mean35, acor_mean35, mean_var35, acov_var35, acor_var35, mean_acov_cor35, mean_acor_cor35, acov_acor_cor35)
est36 = (mean_mean36, acov_mean36, acor_mean36, mean_var36, acov_var36, acor_var36, mean_acov_cor36, mean_acor_cor36, acov_acor_cor36)
est37 = (mean_mean37, acov_mean37, acor_mean37, mean_var37, acov_var37, acor_var37, mean_acov_cor37, mean_acor_cor37, acov_acor_cor37)
est38 = (mean_mean38, acov_mean38, acor_mean38, mean_var38, acov_var38, acor_var38, mean_acov_cor38, mean_acor_cor38, acov_acor_cor38)
est39 = (mean_mean39, acov_mean39, acor_mean39, mean_var39, acov_var39, acor_var39, mean_acov_cor39, mean_acor_cor39, acov_acor_cor39)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
return(tojest)
}
function m_tojmoment(data, acov_order, acor_order, B){
N = rows(data)
S = cols(data)
mean_est = J(N, 1, 0)
acov_est = J(N, 1, 0)
acor_est = J(N, 1, 0)
for (i=1; i<=N; i++){
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
}
if (mod(S,6) == 0) {
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::(S / 3)]
data32 = data[., (S / 3 + 1)::(2*S / 3)]
data33 = data[., (2 * S / 3 + 1)::S]
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
}
equantity = (mean_est, mean_est21, mean_est22, mean_est31, mean_est32, ///
mean_est33, acov_est, acov_est21, acov_est22, acov_est31, acov_est32, ///
acov_est33, acor_est, acor_est21, acor_est22, acor_est31, acor_est32, acor_est33)
estimate_value = tojmomentest0(equantity, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,] = tojmomentest0(equantity, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
} else if (mod(S,6)==1){
// split panel data for T equivalent to 1 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3))]
data33 = data[., (2 * floor(S / 3) + 1)::S]
data34 = data[., 1::floor(S / 3)]
data35 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * floor(S / 3) + 2)::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data39 = data[., (2 * floor(S / 3) + 2)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
equantity = (mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, ///
mean_est36, mean_est37, mean_est38, mean_est39, acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, ///
acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39, acor_est, acor_est21, ///
acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, ///
acor_est38, acor_est39)
estimate_value = tojmomentest1(equantity, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,] = tojmomentest1(equantity, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
} else if (mod(S,6)==2){
// split panel data for T equivalent to 2 modulo 6
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1) ]
data33 = data[., (2 * ceil(S / 3))::S]
data34 = data[., 1::ceil(S / 3)]
data35 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * ceil(S / 3))::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * ceil(S / 3))]
data39 = data[., (2 * ceil(S / 3) + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
equantity = (mean_est, mean_est21, mean_est22, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, ///
mean_est36, mean_est37, mean_est38, mean_est39, acov_est, acov_est21, acov_est22, acov_est31, ///
acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39, acor_est, acor_est21, ///
acor_est22, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, ///
acor_est38, acor_est39)
estimate_value = tojmomentest2(equantity, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,] = tojmomentest2(equantity, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
} else if (mod(S,6) == 3) {
// split panel data for T equivalent to 3 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::(S / 3)]
data32 = data[., (S / 3 + 1)::(2*S / 3)]
data33 = data[., (2 * S / 3 + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
}
equantity = (mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33, ///
acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, ///
acov_est32, acov_est33, acor_est, acor_est21, ///
acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33)
estimate_value = tojmomentest3(equantity, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,] = tojmomentest3(equantity, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
} else if (mod(S,6)==4) {
// split panel data for T equivalent to 4 modulo 6
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3))]
data33 = data[., (2 * floor(S / 3) + 1)::S]
data34 = data[., 1::floor(S / 3)]
data35 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * floor(S / 3) + 2)::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data39 = data[., (2 * floor(S / 3) + 2)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
equantity = (mean_est, mean_est21, mean_est22, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, ///
mean_est36, mean_est37, mean_est38, mean_est39, acov_est, acov_est21, acov_est22, acov_est31, ///
acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39, acor_est, acor_est21, ///
acor_est22, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, ///
acor_est38, acor_est39)
estimate_value = tojmomentest4(equantity, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,] = tojmomentest4(equantity, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
} else {
// split panel data for T equivalent to 5 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1) ]
data33 = data[., (2 * ceil(S / 3))::S]
data34 = data[., 1::ceil(S / 3)]
data35 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * ceil(S / 3))::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * ceil(S / 3))]
data39 = data[., (2 * ceil(S / 3) + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
equantity = (mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, ///
mean_est36, mean_est37, mean_est38, mean_est39, acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, ///
acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39, acor_est, acor_est21, ///
acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, ///
acor_est38, acor_est39)
estimate_value = tojmomentest5(equantity, 1::N)
number_par = length(estimate_value)
estimate_boot = J(B, number_par, 0)
for (b = 1; b <= B; b++) {
index_boot = rdiscrete(N, 1, J(N, 1, 1/N))
estimate_boot[b,] = tojmomentest5(equantity, index_boot)
}
se = sqrt(diagonal(variance(estimate_boot)))
quantile_boot_1 = mm_quantile(estimate_boot, 1, 0.975)
quantile_boot_2 = mm_quantile(estimate_boot, 1, 0.025)
ci_1 = 2 * estimate_value - quantile_boot_1
ci_2 = 2 * estimate_value - quantile_boot_2
ci = (ci_1 \ ci_2)
result = (estimate_value', se, ci')
}
st_matrix("est", estimate_value')
st_matrix("se", se)
st_matrix("ci", ci')
printf("\n")
printf("Estimates for Moments.\n")
printf("Parameters Estimate \n")
printf("________________________________________________________________ \n")
printf("Mean of Mean %f\n",estimate_value[1])
printf("Mean of Autocovariance %f\n",estimate_value[2])
printf("Mean of Autocorrelation %f\n",estimate_value[3])
printf("Variance of Mean %f\n",estimate_value[4])
printf("Variance of Autocovariance %f\n",estimate_value[5])
printf("Variance of Autocorrelation %f\n",estimate_value[6])
printf("Correlation between Mean and Autocovariance %f\n",estimate_value[7])
printf("Correlation between Mean and Autocorrelation %f\n",estimate_value[8])
printf("Correlation between Autocovariance and Autocorelation %f\n",estimate_value[9])
printf("\n")
printf("%f %% Confidence Intervals for Moments.\n", 95)
printf("Parameters Low High\n")
printf("__________________________________________________________________________________\n")
printf("Mean of Mean %f %f\n",ci_1[1],ci_2[1])
printf("Mean of Autocovariance %f %f\n",ci_1[2],ci_2[2])
printf("Mean of Autocorrelation %f %f\n",ci_1[3],ci_2[3])
printf("Variance of Mean %f %f\n",ci_1[4],ci_2[4])
printf("Variance of Autocovariance %f %f\n",ci_1[5],ci_2[5])
printf("Variance of Autocorrelation %f %f\n",ci_1[6],ci_2[6])
printf("Correlation between Mean and Autocovariance %f %f\n",ci_1[7],ci_2[7])
printf("Correlation between Mean and Autocorrelation %f %f\n",ci_1[8],ci_2[8])
printf("Correlation between Autocovariance and Autocorelation %f %f\n",ci_1[9],ci_2[9])
printf("\n")
printf("Standard Errors for Moments.\n")
printf("Parameters Stanadard Errors \n")
printf("________________________________________________________________________\n")
printf("Mean of Mean %f\n",se[1])
printf("Mean of Autocovariance %f\n",se[2])
printf("Mean of Autocorrelation %f\n",se[3])
printf("Variance of Mean %f\n",se[4])
printf("Variance of Autocovariance %f\n",se[5])
printf("Variance of Autocorrelation %f\n",se[6])
printf("Correlation between Mean and Autocovariance %f\n",se[7])
printf("Correlation between Mean and Autocorrelation %f\n",se[8])
printf("Correlation between Autocovariance and Autocorelation %f\n",se[9])
}
/// 4. Kernel Density Estimation
// 4.1. Naive Estimation of Kernel Density
function nekdest (x, X, h) {
N = length(X)
k = normalden((x :- X)/h)
est = sum(k)/(N*h)
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(k:^2) - est^2))
return((est, se))
}
function m_nekd (data, acov_order, acor_order){
N = rows(data)
S = cols(data)
grid = 100
mean_est = J(N,1,0)
acov_est = J(N,1,0)
acor_est = J(N,1,0)
mean_UCI = J(grid,1,0)
acov_UCI = J(grid,1,0)
acor_UCI = J(grid,1,0)
mean_LCI = J(grid,1,0)
acor_LCI = J(grid,1,0)
acov_LCI = J(grid,1,0)
mean_se = J(grid,1,0)
acov_se = J(grid,1,0)
acor_se = J(grid,1,0)
for (i=1 ; i<=N ; i++){
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
}
mean_bw = kdens_bw_dpi(mean_est, level = 2)
acov_bw = kdens_bw_dpi(acov_est, level = 2)
acor_bw = kdens_bw_dpi(acor_est, level = 2)
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_dest = J(grid, 1, 0)
acov_dest = J(grid, 1, 0)
acor_dest = J(grid, 1, 0)
for (i = 1; i <= grid; i++) {
k = nekdest(mean_grid[i], mean_est, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = mean_dest[i] + 1.96 * mean_se[i]
k = nekdest(acov_grid[i], acov_est, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = acov_dest[i] + 1.96 * acov_se[i]
k = nekdest(acor_grid[i], acor_est, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = acor_dest[i] + 1.96 * acor_se[i]
}
temp=st_addvar("double", "mean_dest")
temp=st_addvar("double", "acov_dest")
temp=st_addvar("double", "acor_dest")
temp=st_addvar("double", "mean_grid")
temp=st_addvar("double", "acov_grid")
temp=st_addvar("double", "acor_grid")
temp=st_addvar("double", "mean_LCI")
temp=st_addvar("double", "acov_LCI")
temp=st_addvar("double", "acor_LCI")
temp=st_addvar("double", "mean_UCI")
temp=st_addvar("double", "acov_UCI")
temp=st_addvar("double", "acor_UCI")
st_addobs(max((0,grid - st_nobs())))
st_store(.,"mean_dest", mean_dest\J(st_nobs()-rows(mean_dest),1,.))
st_store(.,"acov_dest", acov_dest\J(st_nobs()-rows(acov_dest),1,.))
st_store(.,"acor_dest", acor_dest\J(st_nobs()-rows(acor_dest),1,.))
st_store(.,"mean_grid", mean_grid\J(st_nobs()-rows(mean_grid),1,.))
st_store(.,"acov_grid", acov_grid\J(st_nobs()-rows(acov_grid),1,.))
st_store(.,"acor_grid", acor_grid\J(st_nobs()-rows(acor_grid),1,.))
st_store(.,"mean_LCI", mean_LCI\J(st_nobs()-rows(mean_LCI),1,.))
st_store(.,"acov_LCI", acov_LCI\J(st_nobs()-rows(acov_LCI),1,.))
st_store(.,"acor_LCI", acor_LCI\J(st_nobs()-rows(acor_LCI),1,.))
st_store(.,"mean_UCI", mean_UCI\J(st_nobs()-rows(mean_UCI),1,.))
st_store(.,"acov_UCI", acov_UCI\J(st_nobs()-rows(acov_UCI),1,.))
st_store(.,"acor_UCI", acor_UCI\J(st_nobs()-rows(acor_UCI),1,.))
}
// 4.2. Half-Panel-Jackknife
function hpjkdest1 (x, X, X1, X2, h) {
N = length(X)
est = sum(normalden((x :- X)/h))/(N * h)
k = normalden((x :- X)/h)
est1 = sum(normalden((x :- X1)/h))/(N * h)
k1 = normalden((x :- X1)/h)
est2= sum(normalden((x :- X2)/h))/(N * h)
k2 = normalden((x :- X2)/h)
hpjest = 2 * est - (est1 + est2) / 2
kest = 2 * k - (k1 + k2) / 2
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(kest:^2) - hpjest^2))
return((hpjest, se))
}
function hpjkdest2 (x, X, X1, X2, X3, X4, h) {
N = length(X)
est = sum(normalden((x :- X)/h))/(N * h)
k = normalden((x :- X)/h)
est1 = sum(normalden((x :- X1)/h))/(N * h)
k1 = normalden((x :- X1)/h)
est2= sum(normalden((x :- X2)/h))/(N * h)
k2 = normalden((x :- X2)/h)
est3 = sum(normalden((x :- X3)/h))/(N * h)
k3 = normalden((x :- X3)/h)
est4= sum(normalden((x :- X4)/h))/(N * h)
k4 = normalden((x :- X4)/h)
hpjest = 2 * est - (est1 + est2 + est3 + est4) / 4
kest = 2 * k - (k1 + k2 + k3 + k4) / 4
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(kest:^2) - hpjest^2))
return((hpjest, se))
}
function m_hpjkd (data, acov_order, acor_order) {
N = rows(data)
S = cols(data)
grid = 100
mean_est = J(N,1,0)
acov_est = J(N,1,0)
acor_est = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
}
mean_bw = kdens_bw_dpi(mean_est, level = 2)
acov_bw = kdens_bw_dpi(acov_est, level = 2)
acor_bw = kdens_bw_dpi(acor_est, level = 2)
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_dest = J(grid, 1, 0)
acov_dest = J(grid, 1, 0)
acor_dest = J(grid, 1, 0)
mean_UCI = J(grid,1,0)
acov_UCI = J(grid,1,0)
acor_UCI = J(grid,1,0)
mean_LCI = J(grid,1,0)
acor_LCI = J(grid,1,0)
acov_LCI = J(grid,1,0)
mean_se = J(grid,1,0)
acov_se = J(grid,1,0)
acor_se = J(grid,1,0)
if (mod(S,2) == 0) {
data1 = data[., 1::(S / 2)]
data2 = data[., (S / 2 + 1)::S]
mean_est1 = J(N,1,0)
acov_est1 = J(N,1,0)
acor_est1 = J(N,1,0)
mean_est2 = J(N,1,0)
acov_est2 = J(N,1,0)
acor_est2 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est1[i] = mean(data1[i,.]')
acov_est1[i] = mataacov(data1[i,.], acov_order)
acor_est1[i] = mataacor(data1[i,.], acor_order)
mean_est2[i] = mean(data2[i,.]')
acov_est2[i] = mataacov(data2[i,.], acov_order)
acor_est2[i] = mataacor(data2[i,.], acor_order)
}
for (i = 1; i <= grid; i++) {
k = hpjkdest1(mean_grid[i], mean_est, mean_est1, mean_est2, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = max((0,mean_dest[i] + 1.96 * mean_se[i]))
mean_dest[i] = max((0,mean_dest[i]))
k = hpjkdest1(acov_grid[i], acov_est, acov_est1, acov_est2, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = max((0,acov_dest[i] + 1.96 * acov_se[i]))
acov_dest[i] = max((0,acov_dest[i]))
k = hpjkdest1(acor_grid[i], acor_est, acor_est1, acor_est2, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = max((0,acor_dest[i] + 1.96 * acor_se[i]))
acor_dest[i] = max((0,acor_dest[i]))
}
}
else {
data1 = data[., 1::floor(S / 2)]
data2 = data[., (floor(S / 2) + 1)::S]
data3 = data[., 1::ceil(S / 2)]
data4 = data[., (ceil(S / 2) + 1)::S]
mean_est1 = J(N,1,0)
acov_est1 = J(N,1,0)
acor_est1 = J(N,1,0)
mean_est2 = J(N,1,0)
acov_est2 = J(N,1,0)
acor_est2 = J(N,1,0)
mean_est3 = J(N,1,0)
acov_est3 = J(N,1,0)
acor_est3 = J(N,1,0)
mean_est4 = J(N,1,0)
acov_est4 = J(N,1,0)
acor_est4 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est1[i] = mean(data1[i,.]')
acov_est1[i] = mataacov(data1[i,.], acov_order)
acor_est1[i] = mataacor(data1[i,.], acor_order)
mean_est2[i] = mean(data2[i,.]')
acov_est2[i] = mataacov(data2[i,.], acov_order)
acor_est2[i] = mataacor(data2[i,.], acor_order)
mean_est3[i] = mean(data3[i,.]')
acov_est3[i] = mataacov(data3[i,.], acov_order)
acor_est3[i] = mataacor(data3[i,.], acor_order)
mean_est4[i] = mean(data4[i,.]')
acov_est4[i] = mataacov(data4[i,.], acov_order)
acor_est4[i] = mataacor(data4[i,.], acor_order)
}
for (i = 1; i <= grid; i++) {
k = hpjkdest2(mean_grid[i], mean_est, mean_est1, mean_est2, mean_est3, mean_est4, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = max((0,mean_dest[i] + 1.96 * mean_se[i]))
mean_dest[i] = max((0,mean_dest[i]))
k = hpjkdest2(acov_grid[i], acov_est, acov_est1, acov_est2, acov_est3, acov_est4, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = max((0,acov_dest[i] + 1.96 * acov_se[i]))
acov_dest[i] = max((0,acov_dest[i]))
k = hpjkdest2(acor_grid[i], acor_est, acor_est1, acor_est2, acor_est3, acor_est4, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = max((0,acor_dest[i] + 1.96 * acor_se[i]))
acor_dest[i] = max((0,acor_dest[i]))
}
}
temp=st_addvar("double", "mean_dest")
temp=st_addvar("double", "acov_dest")
temp=st_addvar("double", "acor_dest")
temp=st_addvar("double", "mean_grid")
temp=st_addvar("double", "acov_grid")
temp=st_addvar("double", "acor_grid")
temp=st_addvar("double", "mean_LCI")
temp=st_addvar("double", "acov_LCI")
temp=st_addvar("double", "acor_LCI")
temp=st_addvar("double", "mean_UCI")
temp=st_addvar("double", "acov_UCI")
temp=st_addvar("double", "acor_UCI")
st_addobs(max((0,grid - st_nobs())))
st_store(.,"mean_dest", mean_dest\J(st_nobs()-rows(mean_dest),1,.))
st_store(.,"acov_dest", acov_dest\J(st_nobs()-rows(acov_dest),1,.))
st_store(.,"acor_dest", acor_dest\J(st_nobs()-rows(acor_dest),1,.))
st_store(.,"mean_grid", mean_grid\J(st_nobs()-rows(mean_grid),1,.))
st_store(.,"acov_grid", acov_grid\J(st_nobs()-rows(acov_grid),1,.))
st_store(.,"acor_grid", acor_grid\J(st_nobs()-rows(acor_grid),1,.))
st_store(.,"mean_LCI", mean_LCI\J(st_nobs()-rows(mean_LCI),1,.))
st_store(.,"acov_LCI", acov_LCI\J(st_nobs()-rows(acov_LCI),1,.))
st_store(.,"acor_LCI", acor_LCI\J(st_nobs()-rows(acor_LCI),1,.))
st_store(.,"mean_UCI", mean_UCI\J(st_nobs()-rows(mean_UCI),1,.))
st_store(.,"acov_UCI", acov_UCI\J(st_nobs()-rows(acov_UCI),1,.))
st_store(.,"acor_UCI", acor_UCI\J(st_nobs()-rows(acor_UCI),1,.))
}
// 4.3. Third-Order-Jackknife
/*computing TOJ kernel density estimate for T equivalent to 0 modulo 6
x point at which the density is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ T/2
X22 vector of half-panel cross-sectional data based on time series (T/2 + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ T/3
X32 vector of one-third-panel cross-sectional data based on time series (T/3 + 1) ~ 2 * T/3
X33 vector of one-third-panel cross-sectional data based on time series 2 * T/3 + 1 ~ T
h bandwidth
*/
function tojkdest0(x, X, X21, X22, X31, X32, X33, h) {
// sample size
N = length(X)
// estimates
est = sum(normalden( (x :- X) / h)) / (N * h)
est21 = sum(normalden( (x :- X21) / h)) / (N * h)
est22 = sum(normalden( (x :- X22) / h)) / (N * h)
est31 = sum(normalden( (x :- X31) / h)) / (N * h)
est32 = sum(normalden( (x :- X32) / h)) / (N * h)
est33 = sum(normalden( (x :- X33) / h)) / (N * h)
k = normalden((x :- X)/h)
k21 = normalden((x :- X21)/h)
k22 = normalden((x :- X22)/h)
k31 = normalden((x :- X31)/h)
k32 = normalden((x :- X32)/h)
k33 = normalden((x :- X33)/h)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33) / 3
kest = 3.536 * k - 4.072 * (k21 + k22) / 2 + 1.536 * (k31 + k32 + k33) / 3
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(kest:^2) - tojest^2))
return((tojest, se))
}
/* computing TOJ kernel density estimate for T equivalent to 1 modulo 6
x point at which the density is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ floor(T/2)
X22 vector of half-panel cross-sectional data based on time series (floor(T/2) + 1) ~ T
X23 vector of half-panel cross-sectional data based on time series 1 ~ ceiling(T/2)
X24 vector of half-panel cross-sectional data based on time series (ceiling(T/2) + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X32 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3))
X33 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 1) ~ T
X34 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X35 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3) + 1)
X36 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 2) ~ T
X37 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X38 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * floor(T/3) + 1)
X39 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 2) ~ T
h bandwidth
*/
function tojkdest1(x, X, X21, X22, X23, X24, X31, X32, X33, X34, X35, X36, X37, X38, X39, h) {
// sample size
N = length(X)
// estimates
est = sum(normalden( (x :- X) / h)) / (N * h)
est21 = sum(normalden( (x :- X21) / h)) / (N * h)
est22 = sum(normalden( (x :- X22) / h)) / (N * h)
est23 = sum(normalden( (x :- X23) / h)) / (N * h)
est24 = sum(normalden( (x :- X24) / h)) / (N * h)
est31 = sum(normalden( (x :- X31) / h)) / (N * h)
est32 = sum(normalden( (x :- X32) / h)) / (N * h)
est33 = sum(normalden( (x :- X33) / h)) / (N * h)
est34 = sum(normalden( (x :- X34) / h)) / (N * h)
est35 = sum(normalden( (x :- X35) / h)) / (N * h)
est36 = sum(normalden( (x :- X36) / h)) / (N * h)
est37 = sum(normalden( (x :- X37) / h)) / (N * h)
est38 = sum(normalden( (x :- X38) / h)) / (N * h)
est39 = sum(normalden( (x :- X39) / h)) / (N * h)
k = normalden((x :- X)/h)
k21 = normalden((x :- X21)/h)
k22 = normalden((x :- X22)/h)
k23 = normalden((x :- X23)/h)
k24 = normalden((x :- X24)/h)
k31 = normalden((x :- X31)/h)
k32 = normalden((x :- X32)/h)
k33 = normalden((x :- X33)/h)
k34 = normalden((x :- X34)/h)
k35 = normalden((x :- X35)/h)
k36 = normalden((x :- X36)/h)
k37 = normalden((x :- X37)/h)
k38 = normalden((x :- X38)/h)
k39 = normalden((x :- X39)/h)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
kest = 3.536 * k - 4.072 * (k21 + k22 + k23 + k24) / 4 + 1.536 * (k31 + k32 + k33 + k34 + k35 + k36 + k37 + k38 + k39) / 9
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(kest:^2) - tojest^2))
return((tojest, se))
}
/* computing TOJ kernel density estimate for T equivalent to 2 modulo 6
x point at which the density is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ T/2
X22 vector of half-panel cross-sectional data based on time series (T/2 + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X32 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3) + 1)
X33 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3)) ~ T
X34 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X35 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * floor(T/3) + 1)
X36 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3)) ~ T
X37 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X38 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * ceiling(T/3))
X39 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3) + 1) ~ T
h bandwidth
*/
function tojkdest2(x, X, X21, X22, X31, X32, X33, X34, X35, X36, X37, X38, X39, h) {
// sample size
N = length(X)
// estimates
est = sum(normalden( (x :- X) / h)) / (N * h)
est21 = sum(normalden( (x :- X21) / h)) / (N * h)
est22 = sum(normalden( (x :- X22) / h)) / (N * h)
est31 = sum(normalden( (x :- X31) / h)) / (N * h)
est32 = sum(normalden( (x :- X32) / h)) / (N * h)
est33 = sum(normalden( (x :- X33) / h)) / (N * h)
est34 = sum(normalden( (x :- X34) / h)) / (N * h)
est35 = sum(normalden( (x :- X35) / h)) / (N * h)
est36 = sum(normalden( (x :- X36) / h)) / (N * h)
est37 = sum(normalden( (x :- X37) / h)) / (N * h)
est38 = sum(normalden( (x :- X38) / h)) / (N * h)
est39 = sum(normalden( (x :- X39) / h)) / (N * h)
k = normalden((x :- X)/h)
k21 = normalden((x :- X21)/h)
k22 = normalden((x :- X22)/h)
k31 = normalden((x :- X31)/h)
k32 = normalden((x :- X32)/h)
k33 = normalden((x :- X33)/h)
k34 = normalden((x :- X34)/h)
k35 = normalden((x :- X35)/h)
k36 = normalden((x :- X36)/h)
k37 = normalden((x :- X37)/h)
k38 = normalden((x :- X38)/h)
k39 = normalden((x :- X39)/h)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
kest = 3.536 * k - 4.072 * (k21 + k22) / 2 + 1.536 * (k31 + k32 + k33 + k34 + k35 + k36 + k37 + k38 + k39) / 9
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(kest:^2) - tojest^2))
return((tojest, se))
}
/* computing TOJ kernel density estimate for T equivalent to 3 modulo 6
x point at which the density is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ floor(T/2)
X22 vector of half-panel cross-sectional data based on time series (floor(T/2) + 1) ~ T
X23 vector of half-panel cross-sectional data based on time series 1 ~ ceiling(T/2)
X24 vector of half-panel cross-sectional data based on time series (ceiling(T/2) + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ T/3
X32 vector of one-third-panel cross-sectional data based on time series (T/3 + 1) ~ 2 * T/3
X33 vector of one-third-panel cross-sectional data based on time series 2 * T/3 + 1 ~ T
h bandwidth
*/
function tojkdest3(x, X, X21, X22, X23, X24, X31, X32, X33, h) {
// sample size
N = length(X)
// estimates
est = sum(normalden( (x :- X) / h)) / (N * h)
est21 = sum(normalden( (x :- X21) / h)) / (N * h)
est22 = sum(normalden( (x :- X22) / h)) / (N * h)
est23 = sum(normalden( (x :- X23) / h)) / (N * h)
est24 = sum(normalden( (x :- X24) / h)) / (N * h)
est31 = sum(normalden( (x :- X31) / h)) / (N * h)
est32 = sum(normalden( (x :- X32) / h)) / (N * h)
est33 = sum(normalden( (x :- X33) / h)) / (N * h)
k = normalden((x :- X)/h)
k21 = normalden((x :- X21)/h)
k22 = normalden((x :- X22)/h)
k23 = normalden((x :- X23)/h)
k24 = normalden((x :- X24)/h)
k31 = normalden((x :- X31)/h)
k32 = normalden((x :- X32)/h)
k33 = normalden((x :- X33)/h)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33) / 3
kest = 3.536 * k - 4.072 * (k21 + k22 + k23 + k24) / 4 + 1.536 * (k31 + k32 + k33) / 3
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(kest:^2) - tojest^2))
return((tojest, se))
}
/* computing TOJ kernel density estimate for T equivalent to 4 modulo 6
x point at which the density is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ T/2
X22 vector of half-panel cross-sectional data based on time series (T/2 + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X32 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3))
X33 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 1) ~ T
X34 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X35 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3) + 1)
X36 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 2) ~ T
X37 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X38 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * floor(T/3) + 1)
X39 vector of one-third-panel cross-sectional data based on time series (2 * floor(T/3) + 2) ~ T
h bandwidth
*/
function tojkdest4(x, X, X21, X22, X31, X32, X33, X34, X35, X36, X37, X38, X39, h) {
// sample size
N = length(X)
// estimates
est = sum(normalden( (x :- X) / h)) / (N * h)
est21 = sum(normalden( (x :- X21) / h)) / (N * h)
est22 = sum(normalden( (x :- X22) / h)) / (N * h)
est31 = sum(normalden( (x :- X31) / h)) / (N * h)
est32 = sum(normalden( (x :- X32) / h)) / (N * h)
est33 = sum(normalden( (x :- X33) / h)) / (N * h)
est34 = sum(normalden( (x :- X34) / h)) / (N * h)
est35 = sum(normalden( (x :- X35) / h)) / (N * h)
est36 = sum(normalden( (x :- X36) / h)) / (N * h)
est37 = sum(normalden( (x :- X37) / h)) / (N * h)
est38 = sum(normalden( (x :- X38) / h)) / (N * h)
est39 = sum(normalden( (x :- X39) / h)) / (N * h)
k = normalden((x :- X)/h)
k21 = normalden((x :- X21)/h)
k22 = normalden((x :- X22)/h)
k31 = normalden((x :- X31)/h)
k32 = normalden((x :- X32)/h)
k33 = normalden((x :- X33)/h)
k34 = normalden((x :- X34)/h)
k35 = normalden((x :- X35)/h)
k36 = normalden((x :- X36)/h)
k37 = normalden((x :- X37)/h)
k38 = normalden((x :- X38)/h)
k39 = normalden((x :- X39)/h)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22) / 2 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
kest = 3.536 * k - 4.072 * (k21 + k22) / 2 + 1.536 * (k31 + k32 + k33 + k34 + k35 + k36 + k37 + k38 + k39) / 9
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(kest:^2) - tojest^2))
return((tojest, se))
}
/* computing TOJ kernel density estimate for T equivalent to 5 modulo 6
x point at which the density is estimated
X vector of original cross-sectional data
X21 vector of half-panel cross-sectional data based on time series 1 ~ floor(T/2)
X22 vector of half-panel cross-sectional data based on time series (floor(T/2) + 1) ~ T
X23 vector of half-panel cross-sectional data based on time series 1 ~ ceiling(T/2)
X24 vector of half-panel cross-sectional data based on time series (ceiling(T/2) + 1) ~ T
X31 vector of one-third-panel cross-sectional data based on time series 1 ~ floor(T/3)
X32 vector of one-third-panel cross-sectional data based on time series (floor(T/3) + 1) ~ (2 * floor(T/3) + 1)
X33 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3)) ~ T
X34 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X35 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * floor(T/3) + 1)
X36 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3)) ~ T
X37 vector of one-third-panel cross-sectional data based on time series 1 ~ ceiling(T/3)
X38 vector of one-third-panel cross-sectional data based on time series (ceiling(T/3) + 1) ~ (2 * ceiling(T/3))
X39 vector of one-third-panel cross-sectional data based on time series (2 * ceiling(T/3) + 1) ~ T
h bandwidth
*/
function tojkdest5(x, X, X21, X22, X23, X24, X31, X32, X33, X34, X35, X36, X37, X38, X39, h) {
// sample size
N = length(X)
// estimates
est = sum(normalden( (x :- X) / h)) / (N * h)
est21 = sum(normalden( (x :- X21) / h)) / (N * h)
est22 = sum(normalden( (x :- X22) / h)) / (N * h)
est23 = sum(normalden( (x :- X23) / h)) / (N * h)
est24 = sum(normalden( (x :- X24) / h)) / (N * h)
est31 = sum(normalden( (x :- X31) / h)) / (N * h)
est32 = sum(normalden( (x :- X32) / h)) / (N * h)
est33 = sum(normalden( (x :- X33) / h)) / (N * h)
est34 = sum(normalden( (x :- X34) / h)) / (N * h)
est35 = sum(normalden( (x :- X35) / h)) / (N * h)
est36 = sum(normalden( (x :- X36) / h)) / (N * h)
est37 = sum(normalden( (x :- X37) / h)) / (N * h)
est38 = sum(normalden( (x :- X38) / h)) / (N * h)
est39 = sum(normalden( (x :- X39) / h)) / (N * h)
k = normalden((x :- X)/h)
k21 = normalden((x :- X21)/h)
k22 = normalden((x :- X22)/h)
k23 = normalden((x :- X23)/h)
k24 = normalden((x :- X24)/h)
k31 = normalden((x :- X31)/h)
k32 = normalden((x :- X32)/h)
k33 = normalden((x :- X33)/h)
k34 = normalden((x :- X34)/h)
k35 = normalden((x :- X35)/h)
k36 = normalden((x :- X36)/h)
k37 = normalden((x :- X37)/h)
k38 = normalden((x :- X38)/h)
k39 = normalden((x :- X39)/h)
// TOJ estimate
tojest = 3.536 * est - 4.072 * (est21 + est22 + est23 + est24) / 4 + 1.536 * (est31 + est32 + est33 + est34 + est35 + est36 + est37 + est38 + est39) / 9
kest = 3.536 * k - 4.072 * (k21 + k22 + k23 + k24) / 4 + 1.536 * (k31 + k32 + k33 + k34 + k35 + k36 + k37 + k38 + k39) / 9
se = sqrt((1 / (N - 1)) * ((1 / (N * h^2)) * sum(kest:^2) - tojest^2))
return((tojest, se))
}
function m_tojkd(data, acov_order, acor_order) {
N = rows(data)
S = cols(data)
grid = 100
mean_est = J(N,1,0)
acov_est = J(N,1,0)
acor_est = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est[i] = mean(data[i,.]')
acov_est[i] = mataacov(data[i,.], acov_order)
acor_est[i] = mataacor(data[i,.], acor_order)
}
mean_bw = kdens_bw_dpi(mean_est, level = 2)
acov_bw = kdens_bw_dpi(acov_est, level = 2)
acor_bw = kdens_bw_dpi(acor_est, level = 2)
mean_lim = (min(mean_est), max(mean_est))
acov_lim = (min(acov_est), max(acov_est))
acor_lim = (min(acor_est), max(acor_est))
mean_grid = rangen(mean_lim[1], mean_lim[2], grid)
acov_grid = rangen(acov_lim[1], acov_lim[2], grid)
acor_grid = rangen(acor_lim[1], acor_lim[2], grid)
mean_dest = J(grid, 1, 0)
acov_dest = J(grid, 1, 0)
acor_dest = J(grid, 1, 0)
mean_UCI = J(grid,1,0)
acov_UCI = J(grid,1,0)
acor_UCI = J(grid,1,0)
mean_LCI = J(grid,1,0)
acor_LCI = J(grid,1,0)
acov_LCI = J(grid,1,0)
mean_se = J(grid,1,0)
acov_se = J(grid,1,0)
acor_se = J(grid,1,0)
if (mod(S,6) == 0) {
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::(S / 3)]
data32 = data[., (S / 3 + 1)::(2*S / 3)]
data33 = data[., (2 * S / 3 + 1)::S]
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
}
for (i = 1; i <= grid; i++) {
k = tojkdest0(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est31, mean_est32, mean_est33, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = max((0,mean_dest[i] + 1.96 * mean_se[i]))
mean_dest[i] = max((0,mean_dest[i]))
k = tojkdest0(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est31, acov_est32, acov_est33, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = max((0,acov_dest[i] + 1.96 * acov_se[i]))
acov_dest[i] = max((0,acov_dest[i]))
k = tojkdest0(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est31, acor_est32, acor_est33, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = max((0,acor_dest[i] + 1.96 * acor_se[i]))
acor_dest[i] = max((0,acor_dest[i]))
}
} else if (mod(S,6)==1){
// split panel data for T equivalent to 1 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3))]
data33 = data[., (2 * floor(S / 3) + 1)::S]
data34 = data[., 1::floor(S / 3)]
data35 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * floor(S / 3) + 2)::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data39 = data[., (2 * floor(S / 3) + 2)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
for (i = 1; i <= grid; i++) {
k = tojkdest1(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, mean_est36, mean_est37, mean_est38, mean_est39, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = max((0,mean_dest[i] + 1.96 * mean_se[i]))
mean_dest[i] = max((0,mean_dest[i]))
k = tojkdest1(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = max((0,acov_dest[i] + 1.96 * acov_se[i]))
acov_dest[i] = max((0,acov_dest[i]))
k = tojkdest1(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, acor_est38, acor_est39, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = max((0,acor_dest[i] + 1.96 * acor_se[i]))
acor_dest[i] = max((0,acor_dest[i]))
}
} else if (mod(S,6)==2){
// split panel data for T equivalent to 2 modulo 6
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1) ]
data33 = data[., (2 * ceil(S / 3))::S]
data34 = data[., 1::ceil(S / 3)]
data35 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * ceil(S / 3))::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * ceil(S / 3))]
data39 = data[., (2 * ceil(S / 3) + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
for (i = 1; i <= grid; i++) {
k = tojkdest2(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, mean_est36, mean_est37, mean_est38, mean_est39, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = max((0,mean_dest[i] + 1.96 * mean_se[i]))
mean_dest[i] = max((0,mean_dest[i]))
k = tojkdest2(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est31, acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = max((0,acov_dest[i] + 1.96 * acov_se[i]))
acov_dest[i] = max((0,acov_dest[i]))
k = tojkdest2(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, acor_est38, acor_est39, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = max((0,acor_dest[i] + 1.96 * acor_se[i]))
acor_dest[i] = max((0,acor_dest[i]))
}
} else if (mod(S,6) == 3) {
// split panel data for T equivalent to 3 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::(S / 3)]
data32 = data[., (S / 3 + 1)::(2*S / 3)]
data33 = data[., (2 * S / 3 + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
}
for (i = 1; i <= grid; i++) {
k = tojkdest3(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = max((0,mean_dest[i] + 1.96 * mean_se[i]))
mean_dest[i] = max((0,mean_dest[i]))
k = tojkdest3(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, acov_est32, acov_est33, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = max((0,acov_dest[i] + 1.96 * acov_se[i]))
acov_dest[i] = max((0,acov_dest[i]))
k = tojkdest3(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = max((0,acor_dest[i] + 1.96 * acor_se[i]))
acor_dest[i] = max((0,acor_dest[i]))
}
} else if (mod(S,6)==4) {
// split panel data for T equivalent to 4 modulo 6
data21 = data[., 1::(S / 2)]
data22 = data[., (S / 2 + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3))]
data33 = data[., (2 * floor(S / 3) + 1)::S]
data34 = data[., 1::floor(S / 3)]
data35 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * floor(S / 3) + 2)::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data39 = data[., (2 * floor(S / 3) + 2)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
for (i = 1; i <= grid; i++) {
k = tojkdest4(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, mean_est36, mean_est37, mean_est38, mean_est39, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = max((0,mean_dest[i] + 1.96 * mean_se[i]))
mean_dest[i] = max((0,mean_dest[i]))
k = tojkdest4(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est31, acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = max((0,acov_dest[i] + 1.96 * acov_se[i]))
acov_dest[i] = max((0,acov_dest[i]))
k = tojkdest4(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, acor_est38, acor_est39, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = max((0,acor_dest[i] + 1.96 * acor_se[i]))
acor_dest[i] = max((0,acor_dest[i]))
}
} else {
// split panel data for T equivalent to 5 modulo 6
data21 = data[., 1::floor(S / 2)]
data22 = data[., (floor(S / 2) + 1)::S]
data23 = data[., 1::ceil(S / 2)]
data24 = data[., (ceil(S / 2) + 1)::S]
data31 = data[., 1::floor(S / 3)]
data32 = data[., (floor(S / 3) + 1)::(2 * floor(S / 3) + 1) ]
data33 = data[., (2 * ceil(S / 3))::S]
data34 = data[., 1::ceil(S / 3)]
data35 = data[., (ceil(S / 3) + 1)::(2 * floor(S / 3) + 1)]
data36 = data[., (2 * ceil(S / 3))::S]
data37 = data[., 1::ceil(S / 3)]
data38 = data[., (ceil(S / 3) + 1)::(2 * ceil(S / 3))]
data39 = data[., (2 * ceil(S / 3) + 1)::S]
// estimated quantities for split panel data
mean_est21 = J(N,1,0)
mean_est22 = J(N,1,0)
mean_est23 = J(N,1,0)
mean_est24 = J(N,1,0)
mean_est31 = J(N,1,0)
mean_est32 = J(N,1,0)
mean_est33 = J(N,1,0)
mean_est34 = J(N,1,0)
mean_est35 = J(N,1,0)
mean_est36 = J(N,1,0)
mean_est37 = J(N,1,0)
mean_est38 = J(N,1,0)
mean_est39 = J(N,1,0)
acov_est21 = J(N,1,0)
acov_est22 = J(N,1,0)
acov_est23 = J(N,1,0)
acov_est24 = J(N,1,0)
acov_est31 = J(N,1,0)
acov_est32 = J(N,1,0)
acov_est33 = J(N,1,0)
acov_est34 = J(N,1,0)
acov_est35 = J(N,1,0)
acov_est36 = J(N,1,0)
acov_est37 = J(N,1,0)
acov_est38 = J(N,1,0)
acov_est39 = J(N,1,0)
acor_est21 = J(N,1,0)
acor_est22 = J(N,1,0)
acor_est23 = J(N,1,0)
acor_est24 = J(N,1,0)
acor_est31 = J(N,1,0)
acor_est32 = J(N,1,0)
acor_est33 = J(N,1,0)
acor_est34 = J(N,1,0)
acor_est35 = J(N,1,0)
acor_est36 = J(N,1,0)
acor_est37 = J(N,1,0)
acor_est38 = J(N,1,0)
acor_est39 = J(N,1,0)
for (i=1 ; i<=N ; i++){
mean_est21[i] = mean(data21[i,.]')
mean_est22[i] = mean(data22[i,.]')
mean_est23[i] = mean(data23[i,.]')
mean_est24[i] = mean(data24[i,.]')
mean_est31[i] = mean(data31[i,.]')
mean_est32[i] = mean(data32[i,.]')
mean_est33[i] = mean(data33[i,.]')
mean_est34[i] = mean(data34[i,.]')
mean_est35[i] = mean(data35[i,.]')
mean_est36[i] = mean(data36[i,.]')
mean_est37[i] = mean(data37[i,.]')
mean_est38[i] = mean(data38[i,.]')
mean_est39[i] = mean(data39[i,.]')
acov_est21[i] = mataacov(data21[i,.], acov_order)
acov_est22[i] = mataacov(data22[i,.], acov_order)
acov_est23[i] = mataacov(data23[i,.], acov_order)
acov_est24[i] = mataacov(data24[i,.], acov_order)
acov_est31[i] = mataacov(data31[i,.], acov_order)
acov_est32[i] = mataacov(data32[i,.], acov_order)
acov_est33[i] = mataacov(data33[i,.], acov_order)
acov_est34[i] = mataacov(data34[i,.], acov_order)
acov_est35[i] = mataacov(data35[i,.], acov_order)
acov_est36[i] = mataacov(data36[i,.], acov_order)
acov_est37[i] = mataacov(data37[i,.], acov_order)
acov_est38[i] = mataacov(data38[i,.], acov_order)
acov_est39[i] = mataacov(data39[i,.], acov_order)
acor_est21[i] = mataacor(data21[i,.], acor_order)
acor_est22[i] = mataacor(data22[i,.], acor_order)
acor_est23[i] = mataacor(data23[i,.], acor_order)
acor_est24[i] = mataacor(data24[i,.], acor_order)
acor_est31[i] = mataacor(data31[i,.], acor_order)
acor_est32[i] = mataacor(data32[i,.], acor_order)
acor_est33[i] = mataacor(data33[i,.], acor_order)
acor_est34[i] = mataacor(data34[i,.], acor_order)
acor_est35[i] = mataacor(data35[i,.], acor_order)
acor_est36[i] = mataacor(data36[i,.], acor_order)
acor_est37[i] = mataacor(data37[i,.], acor_order)
acor_est38[i] = mataacor(data38[i,.], acor_order)
acor_est39[i] = mataacor(data39[i,.], acor_order)
}
for (i = 1; i <= grid; i++) {
k = tojkdest5(mean_grid[i], mean_est, mean_est21, mean_est22, mean_est23, mean_est24, mean_est31, mean_est32, mean_est33, mean_est34, mean_est35, mean_est36, mean_est37, mean_est38, mean_est39, mean_bw)
mean_dest[i] = k[1]
mean_se[i] = k[2]
mean_LCI[i] = max((0,mean_dest[i] - 1.96 * mean_se[i]))
mean_UCI[i] = max((0,mean_dest[i] + 1.96 * mean_se[i]))
mean_dest[i] = max((0,mean_dest[i]))
k = tojkdest5(acov_grid[i], acov_est, acov_est21, acov_est22, acov_est23, acov_est24, acov_est31, acov_est32, acov_est33, acov_est34, acov_est35, acov_est36, acov_est37, acov_est38, acov_est39, acov_bw)
acov_dest[i] = k[1]
acov_se[i] = k[2]
acov_LCI[i] = max((0,acov_dest[i] - 1.96 * acov_se[i]))
acov_UCI[i] = max((0,acov_dest[i] + 1.96 * acov_se[i]))
acov_dest[i] = max((0,acov_dest[i]))
k = tojkdest5(acor_grid[i], acor_est, acor_est21, acor_est22, acor_est23, acor_est24, acor_est31, acor_est32, acor_est33, acor_est34, acor_est35, acor_est36, acor_est37, acor_est38, acor_est39, acor_bw)
acor_dest[i] = k[1]
acor_se[i] = k[2]
acor_LCI[i] = max((0,acor_dest[i] - 1.96 * acor_se[i]))
acor_UCI[i] = max((0,acor_dest[i] + 1.96 * acor_se[i]))
acor_dest[i] = max((0,acor_dest[i]))
}
}
temp=st_addvar("double", "mean_dest")
temp=st_addvar("double", "acov_dest")
temp=st_addvar("double", "acor_dest")
temp=st_addvar("double", "mean_grid")
temp=st_addvar("double", "acov_grid")
temp=st_addvar("double", "acor_grid")
temp=st_addvar("double", "mean_LCI")
temp=st_addvar("double", "acov_LCI")
temp=st_addvar("double", "acor_LCI")
temp=st_addvar("double", "mean_UCI")
temp=st_addvar("double", "acov_UCI")
temp=st_addvar("double", "acor_UCI")
st_addobs(max((0,grid - st_nobs())))
st_store(.,"mean_dest", mean_dest\J(st_nobs()-rows(mean_dest),1,.))
st_store(.,"acov_dest", acov_dest\J(st_nobs()-rows(acov_dest),1,.))
st_store(.,"acor_dest", acor_dest\J(st_nobs()-rows(acor_dest),1,.))
st_store(.,"mean_grid", mean_grid\J(st_nobs()-rows(mean_grid),1,.))
st_store(.,"acov_grid", acov_grid\J(st_nobs()-rows(acov_grid),1,.))
st_store(.,"acor_grid", acor_grid\J(st_nobs()-rows(acor_grid),1,.))
st_store(.,"mean_LCI", mean_LCI\J(st_nobs()-rows(mean_LCI),1,.))
st_store(.,"acov_LCI", acov_LCI\J(st_nobs()-rows(acov_LCI),1,.))
st_store(.,"acor_LCI", acor_LCI\J(st_nobs()-rows(acor_LCI),1,.))
st_store(.,"mean_UCI", mean_UCI\J(st_nobs()-rows(mean_UCI),1,.))
st_store(.,"acov_UCI", acov_UCI\J(st_nobs()-rows(acov_UCI),1,.))
st_store(.,"acor_UCI", acor_UCI\J(st_nobs()-rows(acor_UCI),1,.))
}
mata mlib create lpanelhetero, dir(PERSONAL) replace
mata mlib add lpanelhetero *()
mata mlib index
end