Summary: The expression you're looking for is:

$$ \frac{1}{4} \left[ (III + IZZ + ZIZ + ZZI) + (XXX - XYY - YXY - YYX)\right] $$

where Pauli string notation like $XYX$ denotes $\sigma_1 \otimes \sigma_2 \otimes \sigma_1$, for example.

To start, we need to write the $n$-qubit GHZ state as an operator, namely

\begin{equation} |\psi\rangle\langle\psi | = \frac{1}{2}\left( |0\rangle \langle 0|^{\otimes n} + |0\rangle \langle 1|^{\otimes n} + |1\rangle \langle 0|^{\otimes n} + |1\rangle \langle 1|^{\otimes n}\right) \tag{1} \end{equation}\begin{equation} |\psi^{(n)}\rangle\langle\psi^{(n)} | = \frac{1}{2}\left( |0\rangle \langle 0|^{\otimes n} + |0\rangle \langle 1|^{\otimes n} + |1\rangle \langle 0|^{\otimes n} + |1\rangle \langle 1|^{\otimes n}\right) \tag{1} \end{equation} where we use the notation \begin{equation} A^{\otimes n} = \underbrace{A \otimes A \otimes \dots \otimes A}_{\text{n times}}\tag{2} \end{equation} Taking the pauli matrices to be $\{I, X, Y, Z\}$, we will break this into its diagonal component and off-diagonal component. Looking at the general $n$-qubit case we can rewrite each term of $|\psi\rangle \langle \psi|$ as one of these: \begin{align} |0\rangle \langle 0|^{\otimes n} &= \frac{1}{2^n}(I + Z)^{\otimes n} \tag{3a-d}\\ |1\rangle \langle 1|^{\otimes n} &= \frac{1}{2^n}(I - Z)^{\otimes n}\\ |0\rangle \langle 1|^{\otimes n} &= \frac{1}{2^n}(X + iY)^{\otimes n}\\ |1\rangle \langle 0|^{\otimes n} &= \frac{1}{2^n}(X - iY)^{\otimes n} \end{align}

The terms on the right each look a bit like $(a + b)^n$ and $(a - b)^n$, which we would be able to expand using the Binomial theorem for $a, b \in \mathbb{R}$. However since we are dealing with non-commuting operators, we will have to settle for something that just looks like the binomial theorem and see how much we can simplify.

Permutations for expressing $(A + B)^{\otimes n}$

Consider a permutation operator $S_\pi$ for some permutation $\pi: \{1, \dots, n\}\rightarrow \{1, \dots, n\}$ which has the effect of rearranging the subsystems of its input according to a given permutation $\pi$. For example if $(\pi(a), \pi(b), \pi(c)) = (b, c, a)$ then the action of $S_\pi$ on a separable operator would be \begin{equation} S_\pi (A \otimes B \otimes C) = B \otimes C \otimes A \tag{4} \end{equation}

Then for two operators $A, B$ with $[A, B] \neq 0$ in general we can write \begin{align} (A + B)^{\otimes n} = A^{\otimes n} + \sum_{\pi} S_\pi\left(A^{\otimes (n-1)} \otimes B\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)}\otimes B^{\otimes 2}\right) + \cdots + B^{\otimes n} \tag{5} \\ (A - B)^{\otimes n} = A^{\otimes n} - \sum_{\pi} S_\pi\left(A^{\otimes (n-1)} \otimes B\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)}\otimes B^{\otimes 2}\right) - \cdots + B^{\otimes n}\tag{6} \end{align}\begin{align} (A + B)^{\otimes n} = A^{\otimes n} + \sum_{\pi} S_\pi\left(A^{\otimes (n-1)} \otimes B\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)}\otimes B^{\otimes 2}\right) + \cdots \tag{5} \\ (A - B)^{\otimes n} = A^{\otimes n} - \sum_{\pi} S_\pi\left(A^{\otimes (n-1)} \otimes B\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)}\otimes B^{\otimes 2}\right) - \cdots \tag{6} \end{align}

With the result thatWe can combine (5)-(6) to find \begin{align} \frac{1}{2} \left((A + B)^{\otimes n} + (A - B)^{\otimes n} \right) &= A^{\otimes n} + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)} \otimes B^{\otimes 2}\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-4)} \otimes B^{\otimes 4}\right) + \cdots \\ &= \sum_{t=0}^{\lfloor n/2\rfloor} \sum_\pi S_\pi\left(A^{\otimes (n-2t)} \otimes B^{\otimes 2t}\right) \tag{7} \end{align}\begin{align} \frac{1}{2} &\left((A + B)^{\otimes n} + (A - B)^{\otimes n} \right) \\&= A^{\otimes n} + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)} \otimes B^{\otimes 2}\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-4)} \otimes B^{\otimes 4}\right) + \cdots \\&= \sum_{t=0}^{\lfloor n/2\rfloor} \sum_\pi S_\pi\left(A^{\otimes (n-2t)} \otimes B^{\otimes 2t}\right) \tag{7} \end{align} where we observe that every term $B^{\otimes k}$ with odd $k$ has been eliminated. This lets us write the following:

\begin{align} |0\rangle \langle 0|^{\otimes n} + |1\rangle \langle 1|^{\otimes n} &= \frac{1}{2^{n}}\left[(I + Z)^{\otimes n} + (I - Z)^{\otimes n}\right] \\&= \frac{1}{2^{n-1}} \sum_{t=0}^{\lfloor n/2\rfloor} \sum_\pi S_\pi\left(I^{\otimes (n-2t)} \otimes Z^{\otimes 2t}\right) \tag{8} \end{align}

and$^1$ \begin{align} |0\rangle \langle 1|^{\otimes n} + |1\rangle \langle 0|^{\otimes n} &= \frac{1}{2^n} \left[(X + iY)^{\otimes n} + (X - iY)^{\otimes n} \right] \\&= \frac{1}{2^{n-1}}\sum_{t=0}^{\lfloor n/2\rfloor} \sum_\pi S_\pi\left(X^{\otimes (n-2t)}\otimes (iY)^{\otimes 2t}\right) \\&= \frac{1}{2^{n-1}}\sum_{t=0}^{\lfloor n/2\rfloor} (-1)^t \sum_\pi S_\pi\left(X^{\otimes (n-2t)}\otimes Y^{\otimes 2t}\right) \tag{9} \end{align}

with the final expression becoming the sum of these two parts: \begin{equation} |\psi\rangle\langle\psi | = \frac{1}{2^{n-1}}\sum_{t=0}^{\lfloor n/2\rfloor} \left[ \sum_\pi S_\pi\left(I^{\otimes (n-2t)} \otimes Z^{\otimes 2t}\right) + (-1)^t \sum_\pi S_\pi\left(X^{\otimes (n-2t)} \otimes Y^{\otimes 2t}\right)\right] \tag{10} \end{equation}\begin{equation} |\psi^{(n)}\rangle\langle\psi^{(n)} | = \frac{1}{2^{n-1}}\sum_{t=0}^{\lfloor n/2\rfloor} \left[ \sum_\pi S_\pi\left(I^{\otimes (n-2t)} \otimes Z^{\otimes 2t}\right) + (-1)^t \sum_\pi S_\pi\left(X^{\otimes (n-2t)} \otimes Y^{\otimes 2t}\right)\right] \tag{10} \end{equation}

Solution for $n=3$

Now we can write out 3-qubit case by just enumerating over Pauli strings with the correct parity. In the sum of (8) we will keep the terms $III$, $IZZ$, $ZIZ$, $ZZI$. In the sum of (9) we will keep the terms $XXX$, $XYY$, $YXY$, $YYX$. Plugging into (10) we get

\begin{align} |\psi^{(3)}\rangle\langle\psi^{(3)} | = \frac{1}{4} \left[ (III + IZZ + ZIZ + ZZI) + (XXX - XYY - YXY - YYX)\right] \tag{11} \end{align}

Direct Fidelity Estimation

Looking at Equations (10) and (11) we can see that as the system size grows, the GHZ state will continue to have support on a tiny subset of the $4^n$ Pauli strings: There are only terms containing $Z$ and $I$ together or $X$ and $Y$ together, and even among these a large chunk are missing for the parity constraints shown above.

One consequence of this is that you can estimate the fidelity of a GHZ on $n$ qubits using only $O(n)$ measurement configurations (Gühne, 2007) compared to something exponential in $n$ that you might expect. This is a special case of the more general idea of Direct Fidelity Estimation (DFE) (Flammia, 2011) which demonstrated how fidelity estimation could be improved based on how many of the $4^n$ Pauli strings were missing from the Pauli representation of your state. Of course, (Huang, 2020) brought this way down to $O(1)$ scaling using classical shadows but this is slightly less intuitive than the basic idea.

$^1$ This recovers Equation (5) in (Gühne, 2007)