Hamilton’s approach is easy.
\begin{align*}
\textbf{Date} \quad \quad \textbf{Equation} \qquad\\
\begin{array}{cc}
\hline
0 \quad & \quad y_0 = \phi y_{-1} + \varepsilon_0 \\
1 \quad & \quad y_1 = \phi y_{0} + \varepsilon_1 \\
2 \quad & \quad y_2 = \phi y_{1} + \varepsilon_2 \\
\vdots \quad & \vdots \\
t \quad & \quad y_1 = \phi y_{t-1} + \varepsilon_t \\
\end{array}
\end{align*}
Now express $y_1$ in terms of $y_0$, $y_2$ in terms of $y_1$ etc
\begin{align*}
y_1 &=\phi(\phi y_{-1} + \varepsilon_0) + \varepsilon_1 = \phi^2 y_{-1} + \phi \varepsilon_0 + \varepsilon_1\\
y_2 &=\phi(\phi^2 y_{-1} + \phi \varepsilon_0 + \varepsilon_1) + \varepsilon_2 = \phi^3 y_{-1} +\phi^2 \varepsilon_0 + \phi \varepsilon_1 + \varepsilon_2
\end{align*}
It therefore follows that,
\begin{align*}
y_t= \phi^{t+1}y_{-1}+ \phi^{t} \varepsilon_0 +\phi^{t-1} \varepsilon_1 +\phi^{t-2} \varepsilon_2 + \dots + \phi \varepsilon_{t-1} + \varepsilon_t
\end{align*}
This is the recursive substitution solution to the first-order difference equation with a reference start point (date) of zero. Although it may look slightly odd, we could re-write the equation to include the reference origin;
\begin{align*}
y_{0+t}= \phi^{t+1}y_{0-1}+ \phi^{t} \varepsilon_{0+0} +\phi^{t-1} \varepsilon_{0+1} +\phi^{t-2} \varepsilon_{0+2} + \dots + \phi \varepsilon_{0+t-1} + \varepsilon_{0+t}
\end{align*}
Note that $y_t$ is given at final time $t$ by a recursive equation of order $t+1$ dynamic multiplier ($\phi$) terms and a disturbance term with a dynamic multiplier of 1. We can generalise the above situation in two ways. We could START at time $y_t$ at which $j=0$ (like at time zero above) and look into the future to concluding point $y_{t + j}$ (like $y_t$ above). Or we could start from END point $y_{t + j}$ and look back into the past to beginning point $y_t$ at which $j=0$. Either way, we could then replace subscript $0$ with $t$ and use $j$ to index time intervals such that
\begin{align*}
y_{t+j}= \phi^{j+1}y_{t-1}+ \phi^{j} \varepsilon_{t+0} +\phi^{j-1} \varepsilon_{t+1} +\phi^{j-2} \varepsilon_{t+2} + \dots + \phi^{j-(j-1)} \varepsilon_{t+(j-1)} + \phi^{(j-j)} \varepsilon_{t+j}
\end{align*}
If we let $\phi^{j+1}=\phi^{k}$, it follows that $\phi^{j}=\phi^{k-1}$ and we can re-write the equation
\begin{align*}
y_t= \phi^{k}y_{t-k}+ \phi^{k-1} \varepsilon_{t-(k-1)} +\phi^{k-2} \varepsilon_{t-(k-2)} +\phi^{k-3} \varepsilon_{t-(k-3)} + \dots + \phi \varepsilon_{t-1} + \varepsilon_t
\end{align*}
Observe that if $j$ is sufficiently large, for $|\phi|$ The autoregressive terms become negligible. We may therefore say
\begin{align*}
y_t = \sum_{j=0}^\infty \phi_1^j \varepsilon_{t-j}
\end{align*}