Let $X_t$ be a continuous-time, finite state, real-valued Markov jump process with generator matrix $G$. We can just assume the process takes values in $\{0,1,2,\ldots,m\}$. I think this works more generally, but this is the case I am concerned with. I discretize the process to create discrete-time Markov chain $\{Y_j^n\}_{j=0,1,\ldots}$ with transition matrix $I+\frac1n G$ (ignoring an $O(1/n^2)$ term to make the notation easier). Then we define a continuous-time process $X_t^n$ which can have jumps at the times $\{\frac1n,\frac2n,\ldots...\}$ by $Y_{t}^n=Y_{\lfloor tn\rfloor}$. So $Y^n_t$ is essentially an approximation to the process $X_t$ which only allows jumps at certain discrete times, but that timescale gets refined as $n$ increases.
We can assume the process $X_t$ is irreducible and that it's discretization $Y_j$ is aperiodic if any of that is important, but I don't think it is. Also, I'll always assume both processes start with the same initial distribution, $q$.
The space of sample paths is the set of cadlag step functions, which, with the Skorohod topology and metric, is a complete and separable metric space. See Billingsley's Convergence of Stochastic Processes for all the fine details, but Theorem 12.6 there (p.136 in 2nd edition) seems to tell me that to show $Y_t^n$ converges weakly to $X_t$, it is sufficient to show convergence for finite-dimensional distributions. I am curious about how to go about arguing this.
My attempt: Let $\{t_1,\ldots,t_k\}$ are some time points, that $P_n$ is the distribution of $Y_t^n$, that $P$ the distribution of $X_t$, and that $\pi_{t_1,\ldots,t_k}$ is the projection operator for those time points. I.e. $\pi_{t_1,\ldots,t_k}x=(x_{t_1},\ldots,x_{t_k})\in\mathbb R^k$ (following Billingsley's notation). So I just need to show that $P_n\pi_{t_1,\ldots,t_k}^{-1}\Rightarrow P\pi_{t_1,\ldots,t_k}^{-1}$, which is convergence in distribution now.
Given that the processes only take values in $S=\{0,1,\ldots,m\}$, then I can just take an arbitrary list of values from that set $s_j\in S, j=1,2,...,k$ and show that $P_n\pi_{t_1,\ldots,t_k}^{-1}(s_1,\ldots,s_k)\to P\pi_{t_1,\ldots,t_k}^{-1}(s_1,\ldots,s_k)$ as $n\to\infty$. I only need to consider such a list of states $(s_j)$ because another other measurable subset of $\mathbb R^k$ can be constructed in a product topology sense using sets of the form $(-\infty,a)\subset\mathbb R$. Of course, I'm brushing aside a lot of details here.
I'll attempt this for $k=1$: $$P_n(Y_{t_1}^n=s_1)=\left(q\left(I+\frac1n G\right)^{\lfloor t_1n\rfloor}\right)_{s_1}\to (qe^{t_1G})_{s_1}=P(X_{t_1}=s_1)$$ as $n\to\infty$, with the notation meaning the $s_1^{\text{th}}$ element of the resulting probability vector.
Now for arbitrary $k$, the calculation will be the same, but just taking the increments of time $t_j-t_{j-1}$ and the fact that the evolution is independent among those, and we'll just have starting states $s_{t_{j-1}}$ and ending states $s_{t_j}$. The convergence of the matrix exponentials will work the same as above and the fact that the matrix exponential products will simplify nicely.
Hence the finite-dimensional distributions converge giving us that $Y_t^n$ weakly converges to $X_t$.
So my questions are: (1) Is the above correct? (2) Is it enough to convince someone proficient with these matters? (3) Are there any other references that will help me build out the argument or make it more solid?