Is there a Bowker-like test to assess the symmetry of squared contingency tables with small expected values and a small number of observations? I ask this question because Bowker test (McNemar extension for $k\times k$ tables with $k>2$) requires a large number of observations for fitting well and also requires (as a rule of thumb) that the difference between symmetric cells be $> 5$.
Currently I'm working with $5 \times 5$ tables (pre-post responses to Likert-scale questionnaire items). I would like to know if there were significant changes in the distribution of students answers and if so, verify if the responses were higher or lower in the Likert-scale in comparison of the first stage of the experiment. 
Here I have a small number of observations (25 students) What kind of statistic test can I use to test table symmetry that meets my requirements?
 A: I think this is going to work. Given a contingency table, the distribution of the off-diagonal elements $n_{ij},\,i \neq j$, conditional on the sum of the complementary off-diagonal cells, $N_{ij} = n_{ij} + n_{ji}$, can be written as the product of $K(K - 1)/2$ binominal random variables, that is
$$ P(\mathbf{n}) = \prod_{i < j} \binom{N_{ij}}{n_{ij}} \pi_{ij}^{n_{ij}} ( 1- \pi_{ij})^{n_{ij}},$$
where $\mathbf{n}$ is a vector with elements $n_{ij}$ and $\pi_{ij} = \mathbb{E}(n_{ij}/N_{ij}|N_{ij})$. Under the null hypothesis of complete symmetry, $\pi_{ij} = \pi_{ji} = 1/2$, and thus the permutation distribution is given by
$$ P_{0}(\mathbf{n}) = \prod_{i < j} \binom{N_{ij}}{n_{ij}} \left( \frac{1}{2}\right)^{N_{ij}}.$$
The exact significance test is performed by evaluating
$$ p = \sum_{k \in K} P_{0}(\mathbf{n}),$$
where $K = \{k : P_{0}(\mathbf{n}) < P(\mathbf{n})\}$.
The following R code should do the job.
bowker <- function(mat) {
  n <- dim(mat)[1]
  P_0 <- 1
  P <- 1
  p <- 0
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      N <- mat[i,j] + mat[j,i]
      pi <- mat[i,j] / N
      P_0 <- P_0 * choose(N, mat[i,j]) * (1/2)^N
      P <- P * choose(N, mat[i,j]) * pi^mat[i,j] * (1 - pi)^mat[i,j]
      if (P > P_0) {p <- p + P_0}
    }
  }
  return(list(P_0 = P_0, P = P, p = p))
}

# asymmetric contingency table
mat_asym <- matrix(seq(1,5), ncol = 5, nrow = 5)
mat_asym
bowker(mat_asym)$p

# symmetric contingency table
mat_sym <- matrix(2, ncol = 5, nrow = 5)
mat_sym
bowker(mat_sym)$p

