diff --git a/R/bw_table.R b/R/bw_table.R
new file mode 100644
index 0000000..f5f92d1
--- /dev/null
+++ b/R/bw_table.R
@@ -0,0 +1,67 @@
+# Script for generating table 4 in Cattaneo, Chandak, Jansson and Ma (2022)
+start.time <- Sys.time()
+#parameters
+n = 1000
+num_sims = 1000
+y = 0.0
+x = 0.0
+#true density function
+true_dens =  stats::dnorm(y, mean = 0, sd=1)
+
+#grid for bandwidth range
+h_grid = c(0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5)*0.6
+
+bw_testing = function(s){
+  print(s)
+  # simulating data
+  x_data = matrix(rnorm(n, mean=0, sd=1))
+  y_data = matrix(rnorm(n, mean=x_data, sd=1))
+
+  bw_star=h
+  #density estimation
+  est =  lpcde::lpcde(y_data=y_data, x_data=x_data, y_grid=y, x=x, bw=bw_star, rbc=TRUE)
+
+  # extracting values of interest
+  f_hat = est$Estimate[,3]
+  f_hat_rbc = est$Estimate[,4]
+  bias = (f_hat - true_dens)
+  bias_rbc = (f_hat_rbc - true_dens)
+  sd = est$Estimate[,5]
+  sd_rbc = est$Estimate[,6]
+  rmse = sqrt(sd^2 +bias^2)
+
+  # coverage probabilities
+  ce = 100*ifelse((abs(bias/sd)<=1.96), 1, 0)
+  ce_rbc = 100*ifelse((abs(bias_rbc/sd_rbc)<=1.96), 1, 0)
+
+  #average width
+  aw_pw = 2*1.96*sd
+  aw_pw_rbc = 2*1.96*sd_rbc
+
+  results = matrix(c(bw_star, bias, sd, rmse, ce, ce_rbc, aw_pw, aw_pw_rbc),
+                   ncol = 8)
+}
+
+table_data = matrix(ncol=8)
+# running simulations
+for (i in 1:length(h_grid)){
+ h = h_grid[i]
+ object = parallel::mclapply(1:num_sims, bw_testing, mc.cores =8)
+  int_avg = matrix(0L, nrow = num_sims, ncol = 8)
+  for (l in 1:num_sims){
+    int_avg[l,] = object[[l]][1,]
+  }
+table_data = rbind(table_data, (colSums(int_avg)/num_sims))
+}
+table_data = table_data[2:nrow(table_data), ]
+table_data[, 2] = abs(table_data[,2])
+colnames(table_data) = c("BW", "bias", "sd", "rmse", "WBC CE", "RBC CE",
+                         "WBC AL", "RBC AL")
+rownames(table_data) = h_grid/0.6
+#print table
+table_data
+print(xtable(table_data), file='bw_table.tex')
+
+end.time <- Sys.time()
+time.taken <- end.time - start.time
+time.taken
diff --git a/R/bw_table.tex b/R/bw_table.tex
new file mode 100644
index 0000000..9a9c19b
--- /dev/null
+++ b/R/bw_table.tex
@@ -0,0 +1,22 @@
+% latex table generated in R 4.1.2 by xtable 1.8-4 package
+% Tue Mar 26 14:26:35 2024
+\begin{table}[ht]
+\centering
+\begin{tabular}{rrrrrrrrr}
+  \hline
+ & BW & bias & sd & rmse & WBC CE & RBC CE & WBC AL & RBC AL \\ 
+  \hline
+0.5 & 0.30 & 0.01 & 0.14 & 0.15 & 100.00 & 100.00 & 0.56 & 1.90 \\ 
+  0.6 & 0.36 & 0.01 & 0.08 & 0.10 & 98.00 & 100.00 & 0.32 & 1.08 \\ 
+  0.7 & 0.42 & 0.01 & 0.05 & 0.06 & 96.00 & 100.00 & 0.20 & 0.67 \\ 
+  0.8 & 0.48 & 0.02 & 0.03 & 0.05 & 89.00 & 99.00 & 0.13 & 0.44 \\ 
+  0.9 & 0.54 & 0.02 & 0.02 & 0.04 & 77.00 & 96.00 & 0.10 & 0.32 \\ 
+  1 & 0.60 & 0.02 & 0.02 & 0.03 & 70.00 & 94.00 & 0.07 & 0.23 \\ 
+  1.1 & 0.66 & 0.03 & 0.01 & 0.04 & 38.00 & 86.00 & 0.05 & 0.18 \\ 
+  1.2 & 0.72 & 0.03 & 0.01 & 0.03 & 37.00 & 79.00 & 0.04 & 0.14 \\ 
+  1.3 & 0.78 & 0.04 & 0.01 & 0.04 & 13.00 & 74.00 & 0.03 & 0.11 \\ 
+  1.4 & 0.84 & 0.04 & 0.01 & 0.04 & 10.00 & 69.00 & 0.03 & 0.09 \\ 
+  1.5 & 0.90 & 0.04 & 0.01 & 0.04 & 4.00 & 76.00 & 0.02 & 0.07 \\ 
+   \hline
+\end{tabular}
+\end{table}
diff --git a/R/lpcde_illustration.R b/R/lpcde_illustration.R
index d29d7a3..408b982 100644
--- a/R/lpcde_illustration.R
+++ b/R/lpcde_illustration.R
@@ -10,18 +10,29 @@ n=1000
 x_data = matrix(rnorm(n, mean=0, sd=1))
 y_data = matrix(rnorm(n, mean=x_data, sd=1))
 y_grid = seq(from=-2, to=2, length.out=10)
-
 # density estimation
 #model1 = lpcde::lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=matrix(c(0, 0), ncol=1), bw=0.5)
-model1 = lpcde::lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=0, bw=0.5)
+model1 = lpcde::lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=0, bw=1)
 summary(model1)
 model_reg = lpcde::lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=0, bw=1, nonneg=TRUE, normalize=TRUE)
 summary(model_reg)
 
+set.seed(42)
+n=1000
+x_data = matrix(rnorm(n, mean=0, sd=1))
+y_data = matrix(rnorm(n, mean=x_data, sd=1))
+y_grid = seq(from=-2, to=2, length.out=10)
+# density estimation
+model1 = lpcde::lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=0)
+model_reg = lpcde::lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=0, nonneg=TRUE, normalize=TRUE)
+#plotting densities
+#standard estimate
 plot(x=y_grid, y=model1$Estimate[,3], type="l", lty=1,
      xlab="", ylab="density", ylim=c(0,0.5))
+#regularized estimate
 lines(x=y_grid, y=model_reg$Estimate[,3], lty=2)
-lines(x=y_grid, y=dnorm(y_grid, mean=mean(x_data), sd=0.8), lty=1,col=2)
+#true density
+lines(x=y_grid, y=dnorm(y_grid, mean = mean(y_data)), lty=1,col=2)
 legend('topright',lwd=1, legend=c('f', 'normalized f', 'true f'),lty = c(1, 2, 1), col=c(1,1,2))
 
 #customizing output