Solution:

We will classify all good functions with parameter $\alpha \in \mathbb{F}_p \setminus \{0\}$. We will use $=$ to denote equality modulo $p$ when appropriate. Let the given statement be $P(x, y)$. Then $P(x, 0)$ gives $f(x) f(0) = 2 f(x)$ for all $x \in \mathbb{Z}$. $f$ cannot be the zero function because it does not have minimal period $2016$. Therefore $f(0) = 2$.

$P(x, 1)$ gives $f(x+1) - f(1) f(x) + \alpha f(x-1) = 0$. This means that $f$ satisfies a recurrence of depth two.

We can interpret $f$ as a function $f: \mathbb{Z} \rightarrow \mathbb{F}_{p^2}$, embedding the codomain $\mathbb{F}_p$ inside of $\mathbb{F}_{p^2}$. Let $\tau, \rho$ be the unique roots of $t^2 - f(1) t + \alpha = 0$ in $\mathbb{F}_{p^2}$, so that $\tau + \rho = f(1)$, $\tau \rho = \alpha$. Since $\alpha$ is nonzero, so are $\tau$ and $\rho$. Notice that $f(0) = \tau^0 + \rho^0$ and $f(1) = \tau^1 + \rho^1$. By induction up and down, we see that $f(n) = \tau^n + \rho^n$ for all $n \in \mathbb{Z}$.

We can also easily check that given any value of $f(1)$ and $\alpha$ with $\alpha \not\equiv 0 (\bmod p)$, we get a unique such function. Thus there are $p(p-1)$ total such functions, one corresponding to each quadratic $t^2 - a t + b$, where $b \neq 0$. If $\tau, \rho$ are the roots of the quadratic, the associated function is $f(n) = \tau^n + \rho^n$.

Now if $f$ is $q$-periodic then $f(n) = f(n+q)$ for all $n \in \mathbb{Z}$, so $\tau^n + \rho^n = \tau^{n+q} + \rho^{n+q}$ for all $n \in \mathbb{Z}$ so
$$\tau^n (1 - \tau^q) = \rho^n (\rho^q - 1)$$
We have two cases.

Case 1. $\tau \neq \rho$.

Then $n=0$ and $n=1$ above give $1 = \tau^q = \rho^q$, and this is furthermore sufficient to be $q$-periodic. Thus the period of $f$ is the minimal $q$ such that $\tau^q = \rho^q = 1$, so the lcm of the orders of $\tau, \rho$ as elements of $\mathbb{F}_{p^2}$.

Case 2. $\tau = \rho$.

Then we see that $\tau = \rho = \frac{f(1)}{2}$, so they are in $\mathbb{F}_p$. This automatically means that $f(n) = 2 \left(\frac{f(1)}{2}\right)^n$ for all $n \in \mathbb{Z}$, where $f(1) \neq 0$ is forced (as $\tau \rho = \alpha$). This has period dividing $p-1$. More specifically, the period is the order of $\frac{f(1)}{2}$ as an element of $\mathbb{F}_p$.

In the first case, note that the orders of $\tau, \rho$ dividing $p-1$ implies that $\tau, \rho$ are actually in $\mathbb{F}_p$, so that we can restrict to counting recurrences $f(n) = a^n + b^n$ where $a, b$ are nonzero elements of $\mathbb{F}_p$.

The unordered pair $\{a, b\}$ determines $a+b = f(1)$ and $ab = \alpha$, so it determines a unique $f$ by our earlier observation. We need to count the number of such pairs with $\operatorname{lcm}(\operatorname{ord}_p(a), \operatorname{ord}_p(b)) = p-1$.

Recall that for $d \mid p-1$, there are $\phi(d)$ elements of order $d$ in $\mathbb{F}_p^{\times}$. If $a = b$ we need $\operatorname{ord}_p(a) = p-1$, so there are $\phi(p-1)$ functions. If $S$ is the number of ordered pairs $(a, b)$ with $\operatorname{lcm}(\operatorname{ord}_p(a), \operatorname{ord}_p(b)) = p-1$, then we see that the desired count is $\frac{\phi(p-1) + S}{2}$ (as each unordered pair $\{a, b\}$ is counted twice as an ordered pair except for those with $a = b$).

Now let $n = p-1$. We are counting pairs with the lcm equal to $n$. The number of pairs with respective orders $d_1, d_2$ is $\phi(d_1) \phi(d_2)$ hence our desired sum is $\sum_{d_1, d_2 \mid n, \operatorname{lcm}(d_1, d_2) = n} \phi(d_1) \phi(d_2)$. Now notice that $\phi$ is a multiplicative function, and the condition on the sum splits among the prime powers dividing $n$, so we can choose how $d_1, d_2$ behave on the prime factors of $n$ independently. Therefore if we define
$$\chi(n) = \sum_{d_1, d_2 \mid n, \operatorname{lcm}(d_1, d_2) = n} \phi(d_1) \phi(d_2)$$
then $\chi$ is multiplicative. Finally, if $q$ is prime and $e \geq 1$ we have
$$\begin{aligned} \chi(q^e) & = \sum_{d_1, d_2 \mid q^e, \operatorname{lcm}(d_1, d_2) = q^e} \phi(d_1) \phi(d_2) \\ & = 2 \phi(q^e) (\phi(1) + \phi(q) + \cdots + \phi(q^{e-1})) + \phi(q^e)^2 \\ & = 2 q^{e-1}(q-1) q^{e-1} + q^{2e-2}(q-1)^2 \\ & = q^{2e-2}(q^2 - 1) = q^{2e} \left(1 - \frac{1}{q^2}\right) \end{aligned}$$
This implies that $\chi(n) = n^2 \prod_{q \mid n} \left(1 - \frac{1}{q^2}\right)$ by multiplicativity. Finally, our original answer is
$$\frac{\phi(p-1) + (p-1)^2 \prod_{q \mid p-1} \left(1 - \frac{1}{q^2}\right)}{2}$$
which we now compute for $p = 2017$. We have $p-1 = 2016 = 2^5 \cdot 3^2 \cdot 7$. Thus $\phi(p-1) = 2^4 \cdot 2 \cdot 3 \cdot 6 = 576$ and $\prod_{q \mid 2016} \left(1 - \frac{1}{q^2}\right) = \frac{3}{4} \cdot \frac{8}{9} \cdot \frac{48}{49} = \frac{32}{49}$. Thus the number of functions is
$$\frac{576 + \frac{32}{49} \cdot 2016^2}{2} = \frac{576 + 2^{15} \cdot 3^4}{2} = 288 + 2^{14} \cdot 3^4 = 1327392$$