This post is **almost** an answer to the question. It provides a
method to directly calculate $F(z;t)$ and bypasses the obstacles 
related to computing the expectation values $\langle E_k \rangle_n$.

It's not hard to see that

\begin{equation}
\begin{array}{lll}
\displaystyle H_t(1v) 
&\displaystyle = \, \big(n^2 - t \big) H_t(v)
&\text{if $|v|=n-1$} \\
\displaystyle H_t(2v) 
&\displaystyle = \, \big( 1 - t \big) \big( n^2 - t \big)  H_t(v)
&\text{if $|v|=n-2$} 
\end{array}
\end{equation}

and consequently

\begin{equation}
\langle H_t \rangle_n \ = \
\left\{ 
\begin{array}{ç}
\displaystyle \ \ \ \, {1 \over n} \big( n^2 - t \big)
\langle H_t \rangle_{n-1} \ \ + \\
\displaystyle {n-1 \over n} \big( 1 - t\big) \big(n^2 - t\big)
\langle H_t \rangle_{n-2}
\end{array} \right.
\end{equation}
 
or, after setting 
$\Bbb{f}_t(n) := {1 \over {n!}} \langle H_t \rangle_n$


\begin{equation}
n^2 \Bbb{f}_t(n) \ = \
\left\{
\begin{array}{c} 
\displaystyle \big(n^2 - t \big)\Bbb{f}_t(n-1) \\
+ \\
\displaystyle \big(1 - t \big)\big(n^2 - t \big) \Bbb{f}_t(n-2)
\end{array} \right.
\end{equation}

Of course $F(z;t) = \sum_{n \geq 0} \, \Bbb{f}_t(n) z^n$ and it will satisfy the following second order homogeneous ODE in light of the linear recurrence:

\begin{equation}
\begin{array}{c}
(*) \ \ \displaystyle A(z) {d^2 \over {dz^2}} F(z;t) \ + \
B(z) {d \over {dz}} F(z;t) \ + \
C(z) F(z;t) \ = \ 0 \\
\text{where} \\
\begin{array}{l}
A(z) \ = \ (1-t)z^3 + z^2 - z \\
B(z) \ = \ 5(1-t)z^2 + 3z - 1 \\
C(z) \ = \ (1-t)(4-t)z + (1-t)
\end{array}
\end{array}
\end{equation} 

I don't know how to solve this ODE for general values
of the parameter $t$. 
Nevertheless, by construction, $F(z;0)$ must be 
the generating 
function for the Fibonacci numbers, i.e. 

\begin{equation} 
\begin{array}{l}
\displaystyle F(z;0) 
&\displaystyle = \ 1 + z + 2z^2 + 3z^3 + 5z^4 + \cdots \\
&\displaystyle = \ {1 \over {(1 - z -z^2)}}
\end{array}
\end{equation}

and one can check directly that $(1 - z -z^2)^{-1}$
is indeed a solution to the ODE when $t=0$. When $t=1$ the ODE's 
general solution is $c_1 + c_2 \int z^{-1} \, (z-1)^{-2} \, dz$ where 
$c_1, c_2$ are constants. However, any fibonacci word $u$ with
$|u| > 0$ must contain a box $\Box \in u$ with $\mathrm{h}(\Box)=1$
and so the statistic $H_1(u)$ must vanish.
Combinatorial realities thus force us to select the constant
solution $F(z;1) = 1$.

We can eliminate the first order term in the $(*)$-ODE
by introducing an appropriate phase in the solution $F(z;t) := e^{\varphi(x;t)} \, K(z;t)$. For $\varphi(z;t):= z^{-{1 \over 2}} \big( (1-t)z^2 + z -1\big)^{-1}$ the $(*)$-ODE becomes

\begin{equation}
\begin{array}{c} \displaystyle
(**) \ \ A(z) {d^2 \over {dz^2}} K(z;t) \ + \ D(z) \, K(z;t) \ = \ 0 \\
\text{where} \\
\begin{array}{l}
\displaystyle A(z) \ = \ (1-t)z^3 + z^2 - z \\
\displaystyle D(z) \ = \ {(1-t)(4-t)z^2 + (1-4t)z - 1 \over {4z}} \\
\end{array}
\end{array}
\end{equation}

I'm not sure if this *regauging* will necessarily help; nevertheless I offer it as a potential avenue to crack the nut, so to speak. 
 

The recurrence above is clearly implemented by 
tridiagonal determinants. Specifically $\Bbb{f}_t(n)$ 
equals the initial $n \times n$ principal minor of the following semi-infinite tridiagonal matrix:


\begin{equation}
\begin{pmatrix}
1-t &  1 - {1 \over 4} t & 0 & 0 & \\
t-1 & 1 - {1 \over 4} t &  1 - {1 \over 9} t  & 0 & \\
0 & t-1 &  1 - {1 \over 9} t &  1 - {1 \over 16}t 
& \\ 0 & 0 & t-1 & 1 - {1 \over 16}t & & \\
& & & & \ddots &
\end{pmatrix}
\end{equation}

So the problem of finding $F(z;t)$ can be viewed as a special
case of the general problem of evaluating generating functions
of the sort

\begin{equation}
\sum_{n \geq 0} \, \det (T_n) \, z^n 
\end{equation}

where $T_n$ is the $n \times n$ leading, principal 
submatrix of a infinite $\Bbb{N} \times \Bbb{N}$
tridiagonal matrix $T$. If it were possible to show
that the associated semi-infinite *Hankel matrix*

\begin{equation}
H := \ \begin{pmatrix}
\Bbb{f}_t(0) & \Bbb{f}_t(1) & \Bbb{f}_t(2) & \Bbb{f}_t(3) & \\
\Bbb{f}_t(1) & \Bbb{f}_t(2) & \Bbb{f}_t(3) & \Bbb{f}_t(4) & \\
\Bbb{f}_t(2) & \Bbb{f}_t(3) & \Bbb{f}_t(4) & \Bbb{f}_t(5) & \\
\Bbb{f}_t(3) & \Bbb{f}_t(4) & \Bbb{f}_t(5) & \Bbb{f}_t(6) & \\
& & & & \ddots 
\end{pmatrix}
\end{equation}

were *positive semi-definite* (i.e. all finite, principal 
minors of $H$ are non-negative) then a result of Alan Sokal
(see theorem 2 of https://arxiv.org/pdf/1804.04498.pdf) would
imply that $F(z;t)$ has an expansion as $J$-type continued fraction 

\begin{equation} 
{\alpha_0 \over {1 - z \gamma_0 \ - \ {\displaystyle z^2 \beta_1 \over {\displaystyle 1 - z \gamma_1 \ - \ {z^2 \beta_2 \over {\displaystyle 1 - z \gamma_2 \ - \ {z^2 \beta_3 \over {\ddots}}}}}}}}
\end{equation}

where $\alpha_0 \geq 0$ is a real number and

\begin{equation}
\begin{array}{ll}
\displaystyle \underline{\beta} 
&\displaystyle = \ \big(\beta_1, \ \beta_2, \ \beta_3, \ \dots \big) 
\ \ \text{with $\beta_k \geq 0$ for all $k \geq 1$} \\
\displaystyle \underline{\gamma} 
&\displaystyle = \ \big(\gamma_0, \ \gamma_1, \ \gamma_2, \ \, \dots \big)
\ \ \text{with $\gamma_k \in \Bbb{R}$ for all $k \geq 0$} 
\end{array}
\end{equation} 


When $t=0$ the Hankel matrix $H$ will consists of Fibonacci 
numbers 

\begin{equation}
H := \ \begin{pmatrix}
1 & 1 & 2 & 3 & \\
1 & 2 & 3 & 5 & \\
2 & 3 & 5 & 8 & \\
3 & 5 & 8 & 13 & \\
& & & & \ddots 
\end{pmatrix}
\end{equation}



which is clearly positive semi-definite. As WimC pointed out to me
here (https://math.stackexchange.com/questions/4055848/what-is-the-j-type-continued-fraction-of-the-generating-function-of-the-fibonacc) 
the $J$-type continued fraction of $F(z;0) = \big(1 - z -z^2 \big)^{-1}$
is itself, i.e.

\begin{equation}
{1 \over {1 - 1 \cdot z  \ - \ {\displaystyle 1 \cdot z^2  \over {\displaystyle 1 - 0 \cdot z \ - \ {0 \cdot z^2 \over {\displaystyle 1 - 0 \cdot z \ - \ {0 \cdot z^2 \over {\ddots}}}}}}}}
\end{equation}

Initial calculations using Mathematica endorse the 
claim that $H$ is positive semi-definite for generic values of $t$,
but the initial parameters $\beta_1, \beta_2, \beta_3$
and $\gamma_0, \gamma_1, \gamma_2$ in the $J$-fraction expansion  
do not have a particularly enlightening shape from which to discern
a pattern.

Any help in these matters would be greatly appreciated.

ines.