How can initial conditions be taken into account to calculate a system's terminal value using the final value theorem or some other technique

Ah ha! There's a lot of obfuscation in the problem as stated, plus a bit of over-emphasis on the final value theorem.

$$H(z) = \frac{z^{-1} \phi / (1 - \phi z^{-1})}{1 - z^{-1} \psi / (1 - \psi z^{-1})} \tag 0$$

OK, so far so good.

It arises from the following feedback system

$$Y(z)(1-z^{-1}) = z^{-1}F_{\psi}\cdot Y(z) + z^{-1}F_{\phi}\cdot X(z) \tag{1}$$

where $F_{\alpha}(z)$ is a 1st order high pass filter with parameter $\alpha$ of the form

$$F_{\alpha}(z) = \frac{\alpha (1 - z^{-1})}{1 - \alpha z^{-1}} \tag{2}$$

This was bothering me yesterday, but I didn't notice the essential problem. The general form of the filter described in (2) is first-order, and the delay operators ($z^{-1}$) on the right side of (1) potentially add a second pole.

The original transfer function (0) is in a screwy form, but it appears to be second-order. So the whole problem manages to worm its way through the pole-zero cancellation check by virtue of being stated in such a non-standard way. This question deserves to be an extra-credit problem in a control systems class, or problem 10 out of 10 in a final by a particularly sadistic teacher, or perhaps a routine problem if the instructor believes in homework sets that no one will be able to fully solve.

The problem is that there's a pole-zero cancellation that snuck right by me and everyone else. It's evident in the left-hand side of (1), where the derivative of $y_k$ is the subject of the equation.

So the reason that you can't solve this problem as stated using the final value theorem is that you cannot adequately represent the system using a transfer function. There may be some way to save this within transfer function notation, but I just tried and failed at the first step, so I'm going to do it in state-space.

So, first, shine some light on this obfuscated* problem. I'll do that by substituting (2) into (1):

$$Y(z)(1-z^{-1}) = z^{-1}\frac{\psi(1-z^{-1})}{1 - \psi z^{-1}}Y(z) + z^{-1}\frac{\phi(1-z^{-1})}{1 - \phi z^{-1}}X(z) \tag 3$$

Now the pole-zero cancellation becomes painfully evident: we're expressing a filter of the form $1 - z^{-1}$ on the left-hand side, but when we divide through by $1 - z^{-1}$, it no longer appears in the equation.

Oh joy.

So, drawing this out as a feedback loop (and changing from $z^{-1}$ to $z$ notation, because I'm a control systems guy):

Notation is going to get a bit awkward here. Bear with me. The transfer function in (2) can be realized as a state-space system (where $x$ is a state variable, not the system input):

$$\begin{aligned} x_k = \alpha x_{k-1} + (\alpha - 1)u_k \\ y_k = -x_k + u_k \end{aligned} \tag 4 $$

So (awkwardly) noting that $x_{1,k}$ is the first element of the state vector $\mathbf{x}$ at time $k$, but is not the system input $x_k$, then the 'top' input to the summing junction is (4) with $\alpha = \phi$ and an extra delay:

$$\begin{aligned} \begin{bmatrix}x_{1,k} \\ x_{2,k}\end{bmatrix} = \begin{bmatrix} \phi & 0 \\ 0 & 0\end{bmatrix} \begin{bmatrix}x_{1,k-1} \\ x_{2,k-1}\end{bmatrix} + \begin{bmatrix}\phi - 1 \\ 1\end{bmatrix} x_k \\ u_{1,k} = x_{2, k-1} - x_{1,k-1} \end{aligned} \tag 5 $$

The output of the system, $y_k$, is itself just a state variable (oh joy -- not sarcastically this time). So just let $$y_k = x_{3, k} \tag 6$$.

The filter in the feedback path is again described by (4), with $\alpha = \psi$ and an extra delay:

$$\begin{aligned} x_{4,k} = \psi x_{4,k-1} + (\psi - 1)x_{3, k} \\ u_{2,k} = x_{4,k-1} - x_{3, k-1} \end{aligned} \tag 7 $$

Finally, the integrator in the forward portion of the loop is:

$$x_{2,k} = x_{2,k-1} + u_{1,k} - u_{2,k} \tag 8$$

This goes together into a fourth-order system:

$$\begin{aligned} \mathbf{x}_k = \begin{bmatrix} \phi & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \\ -1 & 1 & 0 & -1 \\ 0 & 0 & \psi - 1 & \psi \end{bmatrix} \mathbf{x}_{k-1} + \begin{bmatrix}\phi - 1 \\ 1 \\ 0 \\ 0\end{bmatrix}x_k \\ y_k = \begin{bmatrix}0 & 0 & 1 & 0\end{bmatrix} \end{aligned} \tag 9 $$

Recast that into $\mathbf{x}_k = A \mathbf{x}_{k-1} + B\,u_k,\,y_k = C\mathbf{x}_k$, and you can test the $A$ and $C$ matrices to see that it is observable, at least for some combinations of $\phi$ and $\psi$ (it is, I checked). Since there's pole-zero cancellation, we can assume that it isn't controllable (you can check that for yourself).

Because it's observable, we can consider the system $$\mathbf{x}_k = A \mathbf{x}_{k-1} + \mathbf{x}_0 \delta(k)$$

The response of this system to any given $\mathbf{x}_0$ is $$Y(z) = C(\mathbf{I}z - A)^{-1}\mathbf{x}_0$$ Now you can find the system response for each initial value in $\mathbf{x}$ using the final value theorem.

* I looked it up -- "obfuscate" literally means "to darken, to put into shadow". So shining a light on the thing really is the antonym of "obfuscate".

After taking some more time to think about it, I think the problem is finally solved. I haven't digested Tim's answer yet, but from what I can see his approach is different. Furthermore, I thought it would be nice to give a closed-form expression for the final value the system's response will reach given values for $x[-1]$ and $y[-1]$.

The final value is given by

$$\bbox[#f8f1ea, 0.6em, border: 0.15em solid #fd8105]{\begin{align}y[\infty]&=\frac{\phi(1-\psi)x[-1]+y[-1]}{(1-\phi)(1-2\psi)},\quad x[-1]\neq 0,\phi\neq 0,\psi\neq 0\\ y[\infty]&=\frac{y[-1]}{1-2\psi},\quad x[-1]=0\textrm{ or }\phi=0\\y[\infty]&=y[-1]+\frac{\phi x[-1]}{1-\phi},\quad\psi=0\end{align}}\tag{1}$$

As pointed out in Tim's answer and in some comments, the problem is the pole-zero cancellation for $z=1$. The transfer function given in the OP correctly describes the system as an LTI system, i.e., with zero initial conditions, but - due to pole-zero cancellation - it doesn't correctly represent the underlying difference equation if initial conditions are non-zero.

Starting from

$$Y(z)(1-z^{-1})=Y(z)z^{-1}\frac{\psi(1-z^{-1})}{1-\psi z^{-1}}+X(z)z^{-1}\frac{\phi(1-z^{-1})}{1-\phi z^{-1}}\tag{2}$$

we can derive the corresponding difference equation without canceling the zero at $z=1$:

$$y[n]=b_1x[n-1]+b_2x[n-2]+b_3x[n-3]-\ldots\\\ldots -a_1y[n-1]-a_2y[n-2]-a_3y[n-3]\tag{3}$$

where the coefficients $b_k$ and $a_k$ depend on the chosen constants $\phi$ and $\psi$.

Now we can transform $(3)$ back to the $\mathcal{Z}$-domain, using the unilateral $\mathcal{Z}$-transform, which takes the initial conditions into account via the rule

$$\begin{align}\mathcal{Z}\big\{y[n-1]\big\}=z^{-1}Y(z)+y[-1]\end{align}\tag{4}$$

This gives an expression for $Y(z)$, the $\mathcal{Z}$-transform of the output signal resulting from a step at the input with initial conditions $x[-1]$ and $y[-1]$.

Now we can use the final value theorem

$$y[\infty]=\lim_{z\to 1}(1-z^{-1})Y(z)\tag{5}$$

which, after expressing the coefficients $a_k$ and $b_k$ in terms of $\phi$ and $\psi$ and after cancelling the common factor $(1-z^{-1})$ in the numerator and denominator of $(1-z^{-1})Y(z)$ results in Eq. $(1)$ for the final value. Note that $(1)$ is only valid if there actually exists a finite final value. This is not always the case because the system can become unstable even though both high-pass filters are stable.

I've run several simulations of the system confirming Eq. $(1)$. I give one example here: $\phi=0.8$, $\psi=0.4$, $x[-1]=1$, $y[-1]=1$. Eq. $(1)$ predicts a final value of $y[\infty]=37$, which is also the result of the simulation: