What you are describing is a (roughly) correct rendition of Monte-Carlo simulation of the true p-value of the one-sample parametric Kolmogorov-Test, and other tests. Although your explanation is roughly correct, you gloss over a few important details. I will give an alternative explanation that I think is a bit clearer, and might aid in assisting you to understand the process in greater detail. I will focus this answer on the parametric Kolmogorov-Smirnov test, but the procedure is essentially the same for other tests that use different test statistics. I will start by giving the general procedure for implementing this kind of test, and then I will answer your specific questions.
The parametric KS test: This test is used when you have and exchangeable sequence of observable values with some unknown underlying distribution, and you want to test whether the underlying distribution is in some specified parametric class $\mathscr{F} = \{ F_\theta | \theta \in \Theta \}$. the null hypothesis for the test is that the observable values are drawn from a distribution in this class, and the alternative hypothesis is that the observable values are drawn from a distribution that is not in this class.
To conduct this test, suppose you observe $n$ data points. You use these data points to form the empirical distribution function $F_n$ and estimate the unknown parameter $\theta$. Regardless of the estimation method you use, your estimator for the unknown parameter will be some function $H: \mathbb{R}^n \rightarrow \Theta$, and the corresponding test statistic for the parametric KS test is:
$$D(\mathbf{x}) = \sup_x| F_n(x)- F_\hat{\theta}(x) |
\quad \quad \quad \quad \quad
\hat{\theta} = H(\mathbf{x}).$$
Substituting the estimator function into the expression for the test statistic gives:
$$D(\mathbf{x}) = \sup_x| F_n(x)- F_{H(\mathbf{x})}(x) |.$$
(Observe that the test statistic $D$ is explicitly a function of the data $\mathbf{x}$, but it is also implicitly dependent on the estimator $H$.) Higher values of this test statistic constitute greater evidence for the alternative hypothesis, so the true p-value of the parametric KS test is:
$$\begin{equation} \begin{aligned}
p(\mathbf{x})
&\equiv \sup_{\theta \in \Theta_0} \mathbb{P}(D(\mathbf{X}) \geqslant D(\mathbf{x}) | \theta) \\[6pt]
&= \sup_{\theta \in \Theta_0} \mathbb{P} \Big( \sup_x| F_n(x)- F_{H(\mathbf{X})}(x) | \geqslant \sup_x| F_n(x)- F_{H(\mathbf{x})}(x) | \Big| \theta \Big). \\[6pt]
\end{aligned} \end{equation}$$
This p-value function is complicated, and even for simple distributional families and estimation methods, it is not usually possible to obtain a simplified expression for the function. This means that it is usual to try to estimate the p-value via Monte-Carlo methods, just as you propose in your question.
Monte-Carlo estimation: As you can see from the above form for the true p-value function, it involves two different kinds of supremums --- one taken over the distributions being compared internally to the probability statement, and one taken over the probability statement for different parameters in the specified class of distributions in our null hypothesis. The explanation in your question refers to calculating the p-value as a sample proportion from Monte-Carlo simulation, but it does not recognise the fact that the true p-value is actually a supremum over a set of probabilities.
In order to simulate effectively for estimation of the p-value in this problem, we simulate $X \sim F_\theta$ indirectly by generating a random variable from an underlying fixed distribution (that does not depend on the parameter) and then we transform this via a deterministic function that depends on the parameter, yielding a random variable with the desired distribution (e.g., by inverse transformation sampling). This allows us to simulate a data vector $\mathbf{X}$ via an underlying random vector with fixed distribution:
$$\mathbf{U} \sim \text{Dist} \quad \quad \quad \quad \quad \mathbf{X} = f(\mathbf{U}, \theta).$$
We can now write our estimated parameter as $\hat{\theta} = G(\mathbf{u},\theta)$ using a compound function $G$, and so the true p-value can now be written as a supremum over a probability pertaining to the random vector $\mathbf{U}$ as:
$$\begin{equation} \begin{aligned}
p(\mathbf{x})
&= \sup_{\theta \in \Theta_0} \mathbb{P} \Big( \sup_x| F_n(x)- F_{G(\mathbf{U}, \theta)}(x) | \geqslant \sup_x| F_n(x)- F_{H(\mathbf{x})}(x) | \Big). \\[6pt]
\end{aligned} \end{equation}$$
The advantage of this step is that we have now simplified the expression for our p-value by ensuring that the probability operation is no longer conditional on $\theta$. Now, for simulations $k=1,...,M$ we generate the random vectors $\mathbf{U}_1,...,\mathbf{U}_M$ from our fixed distribution and we approximate the p-value by:
$$\hat{p}(\mathbf{x}) = \sup_{\theta \in \Theta_0} \frac{1}{M} \sum_{k=1}^M \mathbb{I} \Big( \sup_x| F_n(x)- F_{G(\mathbf{U}_k, \theta)}(x) | \geqslant \sup_x| F_n(x)- F_{H(\mathbf{x})}(x) | \Big). $$
This function is a Monte-Carlo approximation of the true p-value function. It still involves supremums taken over the observable values of the distribution of interest, and the values of the parameter. In practice, these supremums might be replaced either by the outcome of some maximising algorithm (e.g., gradient ascent), or by taking maximums of values over a fine grid of the argument values.
As $M \rightarrow \infty$ this approximating function will converge to the true p-value function. In practice we take $M$ to be some large finite number of simulations, and we can also calculate the approximation over different values of $M$ to see when it is only changing very slowly.
Addressing your specific questions: The above exposition should implicitly answer each of the specific questions you have asked, but I will address them explicitly for completeness.
- Is it possible, for any arbitrary parametric distribution, to properly calculate the p-value for a Kolmogorov-Smirnov test where the parameters of the null distribution are estimated from the data?
If "properly" means getting an exact closed form expression for the p-value, then no. In this type of test the p-value function is extremely complicated (a supremum over a probability of an inequality involving two other suprema!), so we generally estimate the p-value via some kind of Monte-Carlo estimation. This will give you a "proper calculation" for an arbitrary parametric class of distributions, in the sense that the method will converge to the true p-value function as $M \rightarrow \infty$.
- Or does the choice of parametric distribution determine if this can be achieved?
The above method can be implemented for any parametric class of distributions. In practice, the only real difference in having a more complicated parametric class is that the space $\Theta$ may be larger (making it harder to find a supremum over this space) and the corresponding functions $f$ and $G$ used in the simulation may be more complicated, which could present some challenges.
- What about the Anderson-Darling, Cramer von-Mises tests?
Each of these tests uses a different test statistic, so appropriate changes would need to be made to the above functions to reflect this. These would have different forms for the true p-value function and the approximation, but the basic method would essentially be the same.
- What is the general procedure for estimating the correct p-values?
The method presented here is Monte-Carlo estimation, which is what you outlined in your question. There may be other ways to approximate the true p-value function, but this gives a general procedure that can be implemented for an arbitrary class of distributions.
