Closed form for an infinite product
Solution 1:
This formula is an application of the theorem $1.1$ for $n=2$ from your Chamberland & Straub paper (when $a+b=c+d$) :
$$\tag{1}\prod_{k\ge 0}\frac {(k+a)(k+b)}{(k+c)(k+d)}=\frac {\Gamma(c)\Gamma(d)}{\Gamma(a)\Gamma(b)}$$
Observe first that : $\quad\displaystyle n^2-n-1=\left(n-\frac{\sqrt{5}+1}2\right)\left(n+\frac{\sqrt{5}-1}2\right)\;$
An idea is to use $\quad a:=-\dfrac{\sqrt{5}+1}2,\;b:=\dfrac{\sqrt{5}-1}2$.
note that $\;a+b=-1\;$ and take $c=0,\;d=-1$ but we would have $\;0\cdot (-1)\;$ at the denominator so let's shift all these values with an offset $2$ (i.e. set $n:=k+2$) then : $$a:=2-\frac{\sqrt{5}+1}2,\;b:=2+\frac{\sqrt{5}-1}2,\;c:=2,\;d:=1$$ we then get : \begin{align} \prod_{k\ge 0}\frac {(k+2)^2-(k+2)-1}{(k+2)(k+1)}&=\frac {\Gamma(2)\Gamma(1)}{\Gamma\left(2-\frac{\sqrt{5}+1}2\right)\Gamma\left(2+\frac{\sqrt{5}-1}2\right)}\\ &=\frac 1{\left(1-\frac{\sqrt{5}+1}2\right)\left(1+\frac{\sqrt{5}-1}2\right)\Gamma\left(1-\frac{\sqrt{5}+1}2\right)\Gamma\left(1+\frac{\sqrt{5}-1}2\right)}\\ \tag{2}&=-\frac 1{\Gamma\left(1-\frac{\sqrt{5}+1}2\right)\Gamma\left(\frac{\sqrt{5}+1}2\right)} \end{align}
Using the Euler reflexion formula $\displaystyle \Gamma(z)\Gamma(1-z)=\dfrac{\pi}{\sin(\pi z)}$ for $z=\dfrac {\sqrt{5}+1}2$ allows to get : $$\prod_{n\ge 1}\frac {(n+1)^2-(n+1)-1}{n(n+1)}=-\frac {\sin\left(\pi(\sqrt{5}+1)/2\right)}{\pi}$$ and conclude that indeed : $$\tag{3}\prod_{n\ge 1}\left(1-\frac 1{n(n+1)}\right)=-\frac {\sin\left(\pi(\sqrt{5}+1)/2\right)}{\pi}=-\frac1\pi\cos\left(\frac{\sqrt5}{2}\pi\right)$$
Solution 2:
Indeed noticing $n^2-n-1= (n-\phi_+)(n-\phi_-)$ where $\phi_+=\frac{1+\sqrt{5}}{2}$ and $\phi_-=\frac{1-\sqrt{5}}{2}$, We have $$\prod_{n=2}^\infty\left(1-\frac{1}{n(n-1)}\right)=\prod_{n=2}^{\infty}\frac{(n-\phi_+)(n-\phi_-)}{(n-0)(n-1)}=\\ \prod_{n=0}^{\infty}\frac{(n+2-\phi_+)(n+2-\phi_-)}{(n+2)(n+1)}=\frac{\Gamma(2)\Gamma(1)}{\Gamma(2-\phi_+)\Gamma(2-\phi_-)}$$
Where I used the formula $$\prod_{n=0}^{\infty} \frac{(n+a_1)(n+a_2)\cdots(n+a_i)}{(n+b_1)(n+b_2)\cdots(n+b_i)}=\frac{\Gamma(b_1)\Gamma(b_2)\cdots\Gamma(b_i)}{\Gamma(a_1)\Gamma(a_2)\cdots\Gamma(a_i)}$$ which holds whenever $a_1+a_2+\dots+a_i=b_1+b_2+\dots+b_i$.
This follows from Euler's product form of the Gamma function. (and is also the topic of that paper)
The rest follows from the reflection formula $$\Gamma(s)\Gamma(1-s)=\frac{\pi}{\sin(\pi s)}$$, and some adjustments.