Question: Let $X_1,\cdots,X_n \sim \text{IID }\mathcal{N}(\theta,1)$ where $\theta\in\mathbb{R}$ is unknown and let $\psi = \mathbb{P}_\theta(X_1>0)$. Find the maximum likelihood estimator $\hat{\psi}$ of $\psi$ and find a 95% confidence interval for $\psi$.
Comments: I think I understand the first part of the two part question. From the log-likelihood function, it can be shown that $\hat{\theta}_{MLE}=\bar{X}_n=:\sum_{i=1}^n X_i/n.$ Then, with the invariance property of the MLEs, we have $$\hat{\psi}_{MLE} = \mathbb{P}\left(\mathcal{N}(\bar{X}_n,1)>0\right).$$
My question is how do we construct the confidence interval? I am having difficulty calculating the any of the moments $\mathbb{E} \hat{\psi}_{MLE}$ and $\mathbb{E} \hat{\psi}_{MLE}^2$ and hence $\mathbb{V}(\hat{\psi}_{MLE}).$
Since $\bar{X}_n\sim \mathcal{N}(\theta,1/n)$, letting $f_{\bar{X}_n}(x)$ be the density function, I have tried using Tonelli's to switch the order of integration on the following expression $$\mathbb{E} \hat{\psi}_{MLE}=\int_{\mathbb{R}} \mathbb{P}\left(\mathcal{N}(x,1)>0\right) f_{\bar{X}_n}(x)dx=\int_{\mathbb{R}}\int_0^\infty \frac{1}{\sqrt{2\pi}}e^{-(y-x)^2/2}dy \frac{1}{\sqrt{2\pi n}}e^{-(x-\theta)^2/(2\sqrt{n})}dx,$$ but Tonelli's seems to not be the correct route. How do I obtain the variance for the confidence interval?
Edit: Further, I am interested in finding the (nondegenerate) asymptotic distribution of $\hat{\psi}_{MLE}$.