Following Section 7.6 in Diamond & Shurman, the algebraic model of $X(N)$ is constructed (a priori) over $\mathbb{Q}$ by first defining its function field as \begin{equation*} \mathbb{Q}(j,f_{(0,1)}, f_{(1,0)}), \end{equation*} mimicking the complex analytic case. In Exercise 7.7.3 it is claimed that the subgroup which corresponds to this field in the Galois group Gal$(\mathbb{Q}(j,E_j[n])/\mathbb{Q}(j))\simeq$GL$_2(\mathbb{Z}/N \mathbb{Z})$ is just {$\pm I$}, hence explaining why the model is actually defined over $\mathbb{Q}(\mu_n)$.
Why are the matrices \begin{equation*}\pm \begin{bmatrix} 1 & 0 \\ 0 & -1 \\ \end{bmatrix} \end{equation*} not in this subgroup? It seems to me that these also fix both $f_{(1,0)}$ and $f_{(0,1)}$, as they map $P_\tau$ to $\pm P_\tau$ and $Q_\tau$ to $\mp Q_\tau$. Since these matrices have determinant $-1$, the $N$-th roots of unity are not fixed under the Galois action, but instead elements like $\mu_N + \mu_N^{-1}$ are, suggesting that this model of $X(N)$ should be defined over a subfield of $\mathbb{Q}(\mu_N)$ instead.
By letting $R=P+Q$ then the Galois action of one of these matrices on $R$ corresponds to $R^\sigma = P - Q$, suggesting that either $P$ or $Q$ is $2$-torsion, as this action fixes all $x$-coordinates. This is a contradiction for $N>2$, but these matrices still have to define some action on $E_j[N]$ as the full Galois group is GL$_2(\mathbb{Z}/N\mathbb{Z})$, so I am definitely misunderstanding something but I can not see what.
After some thought I have realized why the claim in the last paragraph is wrong. For some $\sigma$ to fix all $x$-coordinates in the torsion subgroup it does not suffice to fix the $x$-coordinates of the generators $P_\tau$ and $Q_\tau$. This is so because the addition formula for points in an elliptic curve is a rational function involving the $x$ and $y$-coordinates of the points to be added, and $\sigma$ need not fix the latter.
Recall that over $\mathbb{C}$ we have that $\mathbb{C}(X(N))$ can be expressed as $\mathbb{C}(j,f_{(0,1)},f_{(1,0)})$ (adjoin $x$-coordinates of the generators of the torsion subgroup) and as $\mathbb{C}(j,\{f_0^{\pm v}\})$ (add the $x$-coordinates of all the torsion points). Since both function fields agree, this suggests that the problem explained in the first paragraph is solved by considering coefficients in $\mathbb{C}$, which in this case means adding $N$-th roots of unity.
That the subgroup in the Galois group corresponding to $\mathbb{Q}(j,f_{(0,1)},f_{(1,0)})$ is $\{\pm I,\pm A\}$, where $\pm A$ are the matrices in the question above, thus seems to be ok, and this just means that $X(N)$ has a model defined over $\mathbb{Q}(\mu_N+\mu_N^{-1})$.
This is also consistent with the model of $X(N)$ defined over $\mathbb{Q}(\mu_N)$. Indeed, consider $\mathbb{Q}(j,\{f_0^{\pm v}\})$. The subgroup in the Galois group corresponding to this function field is $\{\pm I\}$, so it corresponds to a model of $X(N)$ defined over $\mathbb{Q}(\mu_N)$. When we extend scalars to $\mathbb{C}$, both the models of $X(N)$ defined over $\mathbb{Q}(\mu_N+\mu_N^{-1})$ and $\mathbb{Q}(\mu_N)$ recover the function field $\mathbb{C}(X(N))$. Furthermore, when we extend scalars from $\mathbb{Q}(\mu_N+\mu_N^{-1})$ to $\mathbb{Q}(\mu_N)$ the function fields also agree (one can check this is true by degree reasons, but it is also morally true by the reasons explained in the first two paragraphs).