Why does this matrix multiplication converge? [Example of people averaging beliefs] Converging to average
I am pretty new to linear algebra and am working through a problem like this:

*

*There is a true value, $\theta=0$


*Five people have initial guesses for $\theta$ that are drawn from a uniform distribution $[-1, 1]$:


*At each round, people update their guess to an average of the other peoople's guesses


*I found that eventually all people converge to the average of the initial guesses.

Matrix Notation
W is a matrix representing the weights that person $i$ puts on person $j$'s guess.
$G_t$ represents the guesses of each individual at time $t$. $G_t$ = $W^t \times G_0$
$$
W = \begin{pmatrix}
0.00 & 0.25 & 0.25 & 0.25 & 0.25 \\
0.25 & 0.00 & 0.25 & 0.25 & 0.25 \\
0.25 & 0.25 & 0.00 & 0.25 & 0.25 \\
0.25 & 0.25 & 0.25 & 0.00 & 0.25 \\
0.25 & 0.25 & 0.25 & 0.25 & 0.00 \\
\end{pmatrix}
$$
$$ 
G_o = \begin{pmatrix}
\hat{\theta_1}  \\
\hat{\theta_2} \\
\hat{\theta_3}   \\
\hat{\theta_4}  \\
\hat{\theta_5}  \\
\end{pmatrix}
$$
For large $t$, $W^t \times G_0$ approaches the average of the values in $G_0$. Why is that?
Simulation
import numpy as np 
import random 

# True value 
theta = 0

# Number of guessers
n = 5

# Generate guesses around true value
np.random.seed(seed=100)
guesses = np.matrix(np.random.uniform(-1, 1, size=n))

# Weights
mat = np.matrix([[1/(n-1) for x in range(1,n+1)]]*n)
np.fill_diagonal(mat, 0)

# Update beliefs
temp_guess = guesses.copy()
for i in range(10000):
    temp_guess = (np.matmul(mat, temp_guess.T)).T
    
print("Final matrix", temp_guess)
print("Average of initial guesses", np.sum(guesses)/n)

 A: To add to my comment. This only requires Linear Algebra to understand. The eigenvalues in this case can be found in general. Consider $W_n$ where $W^{ij} = 0$ if $i=j$ and $a$ otherwise. Then to find the eigenvalues we have $Wx = \lambda x$, which for each row is
\begin{align}
\sum_{j=1}^n W_{ij}x_j = a \sum_{j=1,j\neq m}^n x_j & = \lambda x_m \\
a \sum_{j=1}^n x_j & = (\lambda+a) x_m.
\end{align}
Since this holds for all $m$ we have $(\lambda+a)x_j = (\lambda + a)x_k$ implying that $\lambda = -a$ or $x_j = x_k$ for all pairs. This implies that another eigenvalue is found by setting $x_j = x$ for all $j$. This leads to $\lambda = a(n-1)$. Finally for your example, note that $a = (n-1)^{-1}$ and so $\lambda = 1$ and $\lambda = -(n-1)^{-1}$. For $W^t$ and $t\to\infty$, the dominant eigenvalue is $\lambda = 1$ with all elements of the eigenvector equal. Hence it will tend towards the average (after normalization).
EDIT:
I suppose it is also important to note that $\lambda = 1$ has multiplicity one. This can be found easily by noting that the trace of the matrix is the sum of the eigenvalues.
SECOND EDIT:
For fun, you can actually generalize the result for all off-diagonal elements to be of size $\beta \in[0,1]$ and main-diagonal elements $\alpha \in [0,1]$ such that $\alpha + (n-1)\beta = 1$ (special case here with $\alpha = 0, \beta = (n-1)^{-1}$). The eigenvalues are $\lambda = 1$ and $\lambda = \alpha - \beta$, with the same convergence to the mean.
THIRD EDIT:
This argument crucially depends on $W = W^T$. If this is not true, then the matrix $W^t$ still converges for $t\to\infty$ but it will tend to the matrix $W_{ij}^t \to Q_{ik} Q_{kj}^{-1}$ where $Q \Lambda Q^{-1} = W$ and $k$ corresponds to the eigenvector with $\lambda = 1$.
