I corrected this post based on the hint below. Results now look consistent.
The derivations follow the steps provided by the appendix of:
“Notes on Estimating the Closed Form of the Hybrid New Phillips Curve”
Jordi Galí, Mark Gertler and J. David López-Salido
Preliminary draft, June 2001
To derive the closed form solution of
x_t = \theta x_{t-1} + (1-\theta)E_t x_{t+1} - \frac{\gamma}{\lambda} \pi_t
define \varepsilon_{t+1} = x_{t+1} -E_t x_{t+1} , substitute E_{t}x_{t+1} by x_{t+1} - \varepsilon_{t+1} and lag the equation by one period
x_{t-1} = \theta x_{t-2} + (1-\theta)(x_{t}-\varepsilon_t) - \frac{\gamma}{\lambda} \pi_{t-1}
Solve for x_t
x_t = \frac{1}{1-\theta} x_{t-1} - \frac{\theta}{1-\theta} x_{t-2} + \frac{\gamma}{\lambda (1-\theta)} \pi_{t-1} + \varepsilon_t
which is a second-order non-homogeneous difference equation in x.
Apply a lag-polynomial
(1- \frac{1}{1-\theta}L + \frac{\theta}{1-\theta}L^2)x_t = \frac{\gamma}{\lambda (1-\theta)} \pi_{t-1} + \varepsilon_t .
Note that 1- \frac{1}{1-\theta}L + \frac{\theta}{1-\theta}L^2 = (1-\delta_1 L)(1-\delta_2 L) holds where \delta_1,\delta_2 are the eigenvalues (see Gandolfo 2009, p.63-64).
Let \mid\delta_2\mid>1 and recall that (1-\delta_2 L) = -\delta_2 L (1-\frac{1}{\delta_2 L}) = -\delta_2 L (1-\frac{1}{\delta_2 }F) where F=L^{-1} is the forward operator.
Thus
\begin{align} (1-\delta_1 L)(1-\delta_2 L)x_t &= \frac{\gamma}{\lambda (1-\theta)} \pi_{t-1} + \varepsilon_t \\
-(1-\delta_1 L)\delta_2 L (1-\frac{1}{\delta_2 }F)x_t &= \frac{\gamma}{\lambda (1-\theta)} \pi_{t-1} + \varepsilon_t\\ (1-\delta_1 L)(1-\frac{1}{\delta_2 }F)x_t &= - \frac{\gamma}{\lambda (1-\theta)\delta_2} \pi_{t} - \frac{1}{\delta_2} \varepsilon_{t+1}
\end{align} .
Recall that (1-\frac{1}{\delta_2 }F)^{-1} = \sum_{i=0}^\infty(\frac{1}{\delta_2})^i F^i.
Thus
\begin{align} (1-\delta_1 L)x_t &= - \frac{\gamma}{\lambda (1-\theta)\delta_2}\sum_{i=0}^\infty\left(\frac{1}{\delta_2}\right)^i \pi_{t+i} -\sum_{i=0}^\infty\left(\frac{1}{\delta_2}\right)^{i+1}\varepsilon_{t+1+i}
\\x_t &= \delta_1 x_{t-1} - \frac{\gamma}{\lambda (1-\theta)\delta_2}\sum_{i=0}^\infty\left(\frac{1}{\delta_2}\right)^i \pi_{t+i} -\sum_{i=0}^\infty\left(\frac{1}{\delta_2}\right)^{i+1}\varepsilon_{t+1+i} \end{align}
Note that the eigenvalues are given by
\begin{align}\delta^2-\frac{1}{1-\theta}\delta + \frac{\theta}{1-\theta} = 0
\\ \delta_{1,2} &= \frac{\frac{1}{1-\theta}\pm\sqrt{(\frac{1}{1-\theta})^2 - 4 \frac{\theta}{1-\theta}}}{2}
\\ &= \frac{1}{1-\theta} \frac{1\pm\sqrt{1 - 4 \theta(1-\theta)}}{2}
\\ \delta_1& = \frac{1}{1-\theta} \frac{1-\sqrt{1 - 4 \theta(1-\theta)}}{2}
\\ &= \frac{\eta}{1-\theta}
\\ \delta_2& = \frac{1}{1-\theta} \frac{1+\sqrt{1 - 4 \theta(1-\theta)}}{2}
\\ &= \frac{1-\eta}{1-\theta} = \frac{1}{\delta_c} \end{align}
where \eta = \frac{1-\sqrt{1 - 4 \theta(1-\theta)}}{2} and \delta_c = \frac{1-\theta}{1-\eta} as in Leitemo 2008, p. 268.
Applying the expectation operator E_t , noting that E_t\varepsilon_{t+1+i} = 0 \forall{} i \geq 0 and plugging in the eigenvalues yields
\begin{align} x_t &= \delta_1 x_{t-1} - \frac{\gamma}{\lambda (1-\theta)\delta_2}\sum_{i=0}^\infty\left(\frac{1}{\delta_2}\right)^i E_t \pi_{t+i}
\\ & = \frac{\eta}{1-\theta} x_{t-1} - \frac{\gamma(1-\theta)}{\lambda (1-\theta)(1-\eta)}\sum_{i=0}^\infty \delta_c^i E_t \pi_{t+i}
\\ & = \frac{\eta}{1-\theta} x_{t-1} - \frac{\gamma}{\lambda (1-\eta)}\sum_{i=0}^\infty \delta_c^i E_t \pi_{t+i}
\\ & = \rho x_{t-1} - \frac{\gamma}{\lambda (1-\eta)}\sum_{i=0}^\infty \delta_c^i E_t \pi_{t+i} \end{align}
which is the closed form solution as in Leitemo (2008), equation (9); since \frac{\eta}{1-\theta} = \frac{\theta}{1-\eta}=\rho according to the result below.