Assume a statistical model $f_{\theta}(X)$ that allows to estimate a parameter vector $\hat{\theta}\in \mathbb{R}^p$ from data $X$ and assume that $p$ is high dimensional (you may assume something like $p=500$).
Since the properties of the estimator cannot be evaluated analytically, I want to show some of them via simulations. Here I draw $N$ times from my data generating process and try to re-estimate the parameters: $\hat{\theta}_{j,i}$ for the $j$th paramter from simulation $i$. My major focus is thereby on unbiasedness, i.e. I want to demonstrate that
$$E[\hat{\theta}]=\theta,$$ This would be a simple task with a low-dimensional parameter vector, where I simply could use boxplots or kernel-density estimates, showing that the distributions of the parameters are centered around the true value. However, this is not possible (or at least extremely confusing) with my high-dimensional parameter vector.
Another issue that complicates matters is given by the fact that the estimation procedure is constructed such that $\sum_{j=1}^p \theta_{j,i} = \sum_{j=1}^p \hat{\theta}_{j,i}$ for $i=1,...,N$ and it is therefore not an option for me to use a measure like $$N^{-1}\sum_{i=1}^N\sum_{j=1}^p( \theta_{j,i} -\hat{\theta}_{j,i}),$$ since it is zero by construction.
Do you have some ideas of how to show unbiasedness in such a setting?
