My answer is in terms of leading questions.
You're trying to establish a fact about what happens with $T_n$ in the limit (as $n\to\infty$). You stopped looking at some finite $n$.
You have demonstrated (more or less) that the gap $|T_n-\theta|$ is small at a specific large $n$ and that it tends to decrease as you approach that $n$.
Did your simulation demonstrate what happens at some larger $n$? Or are you just guessing by extrapolating by eye?
Isn't there a gap between what the simulation suggests and what it actually establishes? To cover that gap, you will end up falling back on some further mathematical argument, will you not? (i.e. "this convergence-suggesting pattern should continue like that at larger $n$ because ...").
Imagine your estimator is a black box, and you only have the simulation to go by.
What makes you sure that the wiggly curve continues to shrink toward $θ$?
Why could it not flatten out some small distance from $θ$ and approach no closer? Why could it not at some later stage increase instead? Has a result at $n=500$ ruled these possibilities out for $n>10^{10^{10}}$?
Ultimately, to decide whether you have satisfied 'consistency', we must look at the definition of consistency -- can a simulation satisfy that definition?
You would need some additional argument (and you may well be able to come up with one) that convinces you that it will indeed continue like you see at larger $n$ -- if such an argument is going to be convincing to a skeptic, it will have elements similar to what is needed for the mathematical argument itself.
That's not to suggest that simulations like this are useless, but there's a difference between a clearly suggestive simulation and a proof of consistency. To span that gap requires some form of argument.
Responding to a comment:
I just want to do a simulation to verify the proof result.
Such a simulation can be useful because if the estimator was inconsistent there might well be a suggestion of that in the results, enough to make you doubt the claim. But not seeing such a suggestion doesn't imply that it's not clear at some much larger $n$.
[Ultimately, consistency may matter less than the properties at sample sizes you might actually observe, but that's not the question we're dealing with here. So let us proceed as if we really do care about asymptotic properties.]
Let's consider some inconsistent estimators. Imagine an estimator that (a) had a bias that didn't shrink toward $0$ as $n$ grew very large; or (b) some limiting non-zero variance.
For (a), consider say estimating a mean, and $T_n=\bar{X}+\delta$ for some tiny $\delta$ (this is obviously biased, but we're pretending we have a black box estimator; imagine it's complicated enough that we can't just see it). At some sample size the variability will have shrunk enough that you'll start to "see" that $|\delta|$ is not zero, so if you're lucky enough to look out that far (and you simulate enough there), you'll have a clear suggestion that this is happening (the estimator will appear to hone in on a value that's not $\theta$).
For (b), consider say $T_n=\frac{6}{\pi^2} (X_1+X_2/4+X_3/9+....)$, which approaches a weighted average of the data (with more weight on the earliest points). At some point - if you simulate out far enough* - you'll start to see that the variance of $T_n$ is not clearly shrinking below some value, at which point you might have cause to hold some doubt as to its consistency.
If you do see effects like these the specific problems that appear can give clues about where it might be productive to focus your attention in looking for an error in the proof.
(You might also like to consider what it might look like when both sorts of effects are present.)
* this particular one doesn't need a large $n$ to see it, but I could make a similar sort of estimator than requires a very large $n$ because the weights shrink slowly enough that you can't clearly see the problem at small $n$
In this Monte-Carlo setup, the indicator function $I: \mathbb{R}^d \rightarrow \{ 0,1 \}$ is defined on the domain $\mathbb{R}^d$, so it already deals with the multivariate nature of the observations and reduces this to a binary output. So, if we have IID observations $\mathbf{x}_1,\mathbf{x}_2,\mathbf{x}_3, ...$ and we set $p_I \equiv \mathbb{P}(I(\mathbf{x}_i) = 1)$ then we have:
$$I(\mathbf{x}_i) \sim \text{IID Bern}(p_I).$$
This then gives the estimator variance:
$$\begin{align}
\mathbb{V}(\hat{\mu})
&= \mathbb{V} \Bigg( \frac{1}{N} \sum_{i=1}^N I(\mathbf{x}_i) \Bigg) \\[6pt]
&= \frac{1}{N^2} \sum_{i=1}^N \mathbb{V} ( I(\mathbf{x}_i) ) \\[6pt]
&= \frac{1}{N^2} \sum_{i=1}^N p_I (1-p_I) \\[6pt]
&= \frac{p_I (1-p_I)}{N}. \\[6pt]
\end{align}$$
As you can see from the working, the reason that this variance does not depend on the dimension $d$ is that $I(\mathbf{x}_i)$ is a Bernoulli random variable irrespective of this dimension and so its variance is determined by $p_I$. So yes, the level of accuracy of the Monte-Carlo estimator is the same whether you have 3 dimensions or 1000. Note that the only thing that you are estimating is the probability of a particular event, which is not something that gets more complicated in higher dimensions.
Best Answer
I'm not sure if it will work 100%, but to show that your estimator has the convergence rate of $O_p(N^{-\frac{1}{3}})$ and you want assess it, from the definition for any given $\epsilon > 0$: $$\lim_{n \rightarrow \infty} Pr(|T_n - \theta| > \epsilon) = 0 $$ So to evaluate this visualize I would suggest fixing some small $\epsilon$ and for each $n$ simulate $|T_n - \theta|$ several times, and calculate the empirical $Pr(|T_n - \theta| > \epsilon) $ and plot it. And with this you can see how the convergence will behave against $N^{-\frac{1}{3}}$. BTW it isn't prove anything, you need to make the theoretical calculations to demonstrate any statement.