Be aware that, when doing MLE (in general, when doing parametric estimation) you are computing (estimating) a parameter of a probability function (pmf). If the variable is discrete, it means (roughly) that its probability function takes discrete values (in this case, $k=1,2,3$), but the parameter itself can be continuous (it can take any real value, in some domain). So, the first thing you need to make clear is that:
what is the parameter of my pmf that I want to estimate? in this case, it's $\theta$
it's continuous? what's its domain? in this case, looking at the pmf, we see that $\theta$ must be in the range $[-1,1]$. In this range, and only in this range the probability function is valid (takes non-negative values). Then the parameter is continous and its domain is $-1 \le\theta \le 1$
Once you have that establlished, you try to write the likelihood. If you are not sure, start by some simple example. Assume you have only two samples, say, $x_1=2$, $x_2=0$. The likelihood of this realization is
$$L(\theta)=p(x_1=2;\theta) \times p(x_2=0;\theta) = \frac{1+\theta}{3} \frac{1-\theta}{3} $$
To write this in general, suppose you have $n_0$ samples that take value $x=0$, $n_1$ that take value $x=1$ etc. Then
$$L(\theta)=p(x=0;\theta)^{n_0}p(x=1;\theta)^{n_1}p(x=2;\theta)^{n_2}$$
Write that expression down, and take its logarithm if you think this simplifies things (it does). Then ask yourself: for given $n_0,n_1,n_2$, this is a (continous) function of $\theta$, what is the value of $\theta$ that maximizes
this function, in the given domain?
Update: given that you've done your homework, here's my solution
$$\log L(\theta)= n_0 \log(1+\theta) +n_2 \log(1-\theta) +\alpha $$
where $\alpha $ is a term that does not depend on $\theta$ (we can leave it out). This function is differentiable in $(-1,1)$, so we can look for critical points (candidate extrema) as:
$$\frac{d\log L(\theta)}{d \theta}= \frac{n_0}{1+\theta}-\frac{n_2}{1-\theta} $$
Equalling this to zero, we get $\theta_0=(n_0-n_2)/(n_0+n_2)$
Have we already found then the MLE? Not really. We have only found a critical point of $L(\theta)$. To assert that a critical point is a global maximum we need to 1) check that it's a local maximum (it could be a local minimum or neither) 2) check that the local maximum is really a global maximum (what about the non-differentiable or boundary points?).
We can usually check that with the second derivative. But in this case it's simpler. We see that at the boundary ($\theta = \pm 1$) the likelihood tends to $-\infty$. Hence, given that the function is differentiable inside the interval, and it has a single critical point, it must be a (local and global) maximum.
I believe I have figured this out, and so am posting it as an answer to my own question.
Using the likelihood function provided in the question above:
$$ \mathcal{L}(\boldsymbol{Z}; \theta) = \left(\frac{1}{2 \pi N \sigma^2}\right)^N \exp\left[-\frac{1}{2 N \sigma^2} \sum_n\left((X_n - p_n)^2 + (Y_n - q_n)^2\right)\right]$$
Since, there are $3$ unknown parameters, $\boldsymbol{\theta} = \begin{bmatrix}t & a & \phi\end{bmatrix}^\textsf{T}$, we must determine the elements of a $3 \times 3$ Fisher Information matrix, $\mathcal{I}$, using:
$$\mathcal{I}_{ij} = -\operatorname{E}\left[
\frac{\partial^2}{\partial\theta_i\, \partial\theta_j} \log \mathcal{L}(\boldsymbol{Z}; \theta)\right] = \frac{1}{N \sigma^2} \sum_n \left[\frac{\partial p_n}{\partial \theta_i} \frac{\partial p_n}{\partial \theta_j} + \frac{\partial q_n}{\partial \theta_i} \frac{\partial q_n}{\partial \theta_j}\right]$$
where $n = -\frac{N}{2}, ..., \frac{N}{2} - 1$.
With this, the $i$th diagonal element of $\mathcal{I}^{-1}$ is the Cramér–Rao lower bound on the variance of an unbiased estimator of $\theta_i$. Using Mathematica to invert $\mathcal{I}$, we get:
$$\operatorname{var}(\hat{t}) = \frac{3 \sigma^2}{\pi^2 a^2} \frac{N^2}{N^2 - 1}$$
$$\operatorname{var}(\hat{a}) = \sigma^2 $$
$$\operatorname{var}(\hat{\phi}) = \frac{\sigma^2}{a^2} \frac{N^2 + 2}{N^2 - 1} $$
I have also experimentally confirmed that the maximum-likelihood estimator outlined in the question seems to achieve these lower bounds.
Best Answer
As I mentioned in a comment, the likelihood can be written in the following convenient way:
\begin{align} L(\theta\mid \boldsymbol x)&=\left(\prod_{i: x_i=0}\frac14\right)\left(\prod_{i:x_i=1}\frac{3\theta}4\right)\left(\prod_{i:x_i=2}\frac{3(1-\theta)}4\right)\quad,\,\theta \in [0,1] \\\\&=\left(\frac14\right)^{n_1}\left(\frac{3\theta}4\right)^{n_2}\left(\frac{3(1-\theta)}4\right)^{n_3}\,, \end{align}
where $n_1,n_2$ and $n_3$ are respectively the number of observations equal to $0,1$ and $2$ in the sample, such that $n_1+n_2+n_3=n$. These $n_i$'s are of course based on the sample observations, which can be explicitly seen by writing them as sums of indicator variables:
$$n_1(\boldsymbol x)=\sum_{i=1}^n \mathbf1(x_i=0)\,,\,n_2(\boldsymbol x)=\sum_{i=1}^n \mathbf1(x_i=1)\,,\,\text{etc.}$$
In fact, it is clear from the likelihood that $(n_1(\boldsymbol X),n_2(\boldsymbol X),n_3(\boldsymbol X))$ has a multinomial distribution.
This should help you in finding the Fisher information.