what you are missing is the automorphism group of the quadratic form $x^2 - 4 x y + y^2.$ On this site, the relevant observation is usually referred to as Vieta Jumping, which is a special case. The point is that a solution $(x,y)$ to $x^2 - 4 x y + y^2= k$ (in this case $k=-2$) can be replaced by a new solution
$$ (y, 4y - x). $$
CHECK!!!!!!!!
This means that $(a_n, a_{n+1})$ becomes $(a_{n+1}, a_{n+2})$ with the specific
$$ a_{n+2} = 4 a_{n+1} - a_n. $$ Thus, they are all integers.
Furthermore, as long as $a_{n+1} > a_n,$ then $a_{n+2} > a_{n+1},$ so they continue increasing and positive forever by induction.
next day: a short tutorial. Suppose we have positive integers $x,y,q$ and a fixed target $T,$ an integer. Then suppose we have
$$ \color{blue}{ x^2 - q xy + y^2 = T.} $$
Let us take $0 < x < y.$ The increasing direction in Vieta Jumping is
$$ \color{green}{ (x,y) \mapsto (y, qy-x).} $$
The decreasing (with $0 < x < y$) direction is
$$ \color{red}{ (x,y) \mapsto ( qx-y, x).} $$
To be specific, the thing keeps decreasing the sum of the entries as long as $qx < 2y.$ The way proofs are arrived at in these problems is to examine those solutions with $x \leq y$ but $qx \geq 2y.$ These are what Hurwitz, in 1907, called Grundlösungen, or fundamental solutions. If there are no fundamental solutions, there are no solutions at all, because positive integers cannot decrease forever. If the fundamental solution has one of the variables $x,y$ equal to zero, that says that the target $T$ is a square. And so on. By and large, the jumping phenomenon stays the same, it is the inequalities towards the end that vary.
Lemma: $f(k)>k$.
Pf: it is clear that $f(1)≥1$ and that $f(1)\neq 1$ so the claim is true for $k=1$. Inductively suppose it is true up to $k-1$. Then $f(k-1)>k-1$ so $f(k-1)≥k$. Since $f(k)>f(k-1)$ we are done.
Claim: $f(n)=n+1$
Pf: Indeed the Lemma shows that $f(n)≥n+1$. But $f(f(n))=n+2\implies f(n)≤n+1$ and we are done.
Best Answer
Let $f:\mathbb{Z}_{>0}\to\mathbb{R}$ be a function that satisfies $$f(xy)+f(x+y)=f(x)\,f(y)+1\tag{*}$$ for all $x,y\in\mathbb{Z}_{>0}$. Let $k:=f(1)$. Plugging in $x:=n$ and $y:=1$ into (*), where $n$ is a positive integer, we get $$f(n+1)=(k-1)\,f(n)+1$$ for every positive integer $n$.
If $k=1$, then $$f(n)=1\text{ for every positive integer $n$}\,.\tag{#}$$ If $k=2$, then we see that $f(n+1)=f(n)+1$ for all positive integers $n$, whence $$f(n)-2=f(n)-f(1)=\sum_{j=1}^{n-1}\,\big(f(j+1)-f(j)\big)=\sum_{j=1}^{n-1}\,1=n-1\,,$$ so that $$f(n)=n+1\text{ for every positive integer $n$}\,.\tag{@}$$ We now assume that $k\notin\{1,2\}$.
Write $g(n):=f(n)+\dfrac{1}{k-2}$ for all $n\in\mathbb{Z}_{>0}$. We then have $$\begin{align}g(n+1)&=f(n+1)+\frac{1}{k-2}=(k-1)\,f(n)+1+\frac{1}{k-2}\\&=(k-1)\,f(n)+\frac{k-1}{k-2}=(k-1)\,g(n)\end{align}$$ for every $n\in\mathbb{Z}_{>0}$. This shows that $$\begin{align}g(n)&=(k-1)^{n-1}\,g(1)=(k-1)^{n-1}\,\left(f(1)+\frac{1}{k-2}\right)\\&=(k-1)^{n-1}\,\left(k+\frac{1}{k-2}\right)=\frac{(k-1)^{n+1}}{k-2}\end{align}$$ for every $n\in\mathbb{Z}_{>0}$. Ergo, $$f(n)=g(n)-\frac{1}{k-2}=\frac{(k-1)^{n+1}-1}{k-2}=t^n+t^{n-1}+\ldots+t+1$$ for every $n\in\mathbb{Z}_{>0}$, where $t:=k-1$.
Now, we have $$f(2)=t^2+t+1$$ and $$f(4)=t^4+t^3+t^2+t+1\,.$$ Substituting $2$ for both $x$ and $y$ in (*) leads to $$2\,f(4)=\big(f(2)\big)^2+1\,.$$ Then, $$2\,\big(t^4+t^3+t^2+t+1\big)=(t^2+t+1)^2+1=t^4+2t^3+3t^2+2t^2+2\,.$$ That is, $t^4=t^2$. Because $k\notin\{1,2\}$, we get $t\notin\{0,1\}$. Hence, $t=-1$ is the only possibility. That is, $k=0$. Consequently, $$f(n)=\frac{1+(-1)^n}{2}=1-(n\,\text{mod}\,2)\text{ for all positive integers $n$}\,.\tag{\$}$$ It is easy to see that the functions (#), (@), and (\$) satisfy the functional equation (*).