Let $g_n = f_n + 1$ so that \begin{align} g_n &= g_{n−1} + g_{n−3} &&\text{for $n\ge3$} \\ g_0 &= 2 \\ g_1 &= 3 \\ g_2 &= 4 \end{align} Let $G(z)=\sum_{n\ge 0} g_n z^n$ be the generating function. The recurrence relation and boundary conditions imply that $$G(z) = 2z^0+3z^1+4z^2+\sum_{n\ge 3} (g_{n−1} + g_{n−3}) z^n = 2+3z+4z^2 + z (G(z)-2-3z) + z^3 G(z),$$ so $$G(z) = \frac{2+z+z^2}{1-z-z^3},$$ which yields $$F(z) = \sum_{n\ge 0} f_n z^n = \sum_{n\ge 0} (g_n-1) z^n = \frac{2+z+z^2}{1-z-z^3} - \frac{1}{1-z} = \frac{1}{(1-z)(1-z-z^3)}.$$

See https://oeis.org/A077868. The linked entry https://oeis.org/A000930 for the shifted sequence $f_{n+3}-1$ has this formula:

a(n) = floor(d*c^n + 1/2) where c is the real root of x^3-x^2-1 and d is the real root of 31*x^3-31*x^2+9*x-1 (c = 1.465571... = A092526 and d = 0.611491991950812...). - Benoit Cloitre, Nov 30 2002