Expectation given pairwise covariances I have 4 variables A,B,C,D over {-1,1} (Rademacher variables) and know that E[A]=E[B]=E[C]=E[D]=0. I also know E[AB], E[BC], E[CD], and E[AD] (but not E[AC], E[BD]), which are all positive. Is there any way to compute E[ABCD]? If not, is there any way to bound it or its magnitude?
 A: This is a linear program.  The variables are the 16 non-negative probabilities for the 16 possible values of $(A,B,C,D)$ subject to linear constraints (sum to unity and the specified expectations).
We might hope the structure of these constraints yields a simple solution, for then it would be worthwhile to derive it analytically.  To this end, we can explore the landscape with the computer.  The expectations must all lie within the range $[-1,1]$ of the values of these variables.  The constraints have a cyclic structure, meaning they have the same form when $(A,B,C,D)$ is cyclically permuted.  We may therefore arrange for $E[AB]$ to be the largest among them.  This reduces the effort to performing a brute-force grid search of the vectors $(E[AB], E[BC], E[CD], E[AD]),$ evaluating the extreme possible values in each case.
It took two minutes to do this for a grid with $1/8$ spacing ($23409$ separate problem specifications).  This scatterplot matrix summarizes the results.  The gray areas plot the values.

It is not bad; but it is sufficiently complicated that I, for one, do not want to bother with further analysis absent some compelling application.  I recommend a one-off solution for your particular case.  Example R code is below.  In the test example the four expected products were set to $1/2, 3/8, 1/4,$ and $1/8.$  The first line of the summary displays all the relevant expectations:
[1] "Min"
   E[A]    E[B]    E[C]    E[D]   E[AB]   E[BC]   E[CD]   E[DA] E[ABCD] 
  0.000   0.000   0.000   0.000   0.500   0.375   0.250   0.125  -0.250 

The last one at right is the minimum value of $E[ABCD].$  Similarly, the maximum is $0.750.$  The rest of the output shows a distribution that realizes the minimum (as a four-tensor) and another that realizes the maximum.
#
# Function to optimize A*B*C*D subject to constraints on all other product
# expectations for variables supported on the set `pm`.
#
f <- function(ab=NA, bc=NA, cd=NA, da=NA, ac=NA, bd=NA,
              abc=NA, abd=NA, acd=NA, bcd=NA,
              a = 0, b = 0, c = 0, d = 0, pm = c(-1, 1)) {
  # The sample space, by rows
  X <- expand.grid(A=pm, B=pm, C=pm, D=pm)
  rownames(X) <- apply(X, 1, paste, collapse = " ")

  # Sum to unity condition
  X <- cbind(Sum = 1, X)

  # Known expectations
  X <- with(X, cbind(X, AB=A*B, BC=B*C, CD=C*D, DA=D*A, AC=A*C, BD=B*D,
                     ABC=A*B*C, ABD=A*B*D, ACD=A*C*D, bcd=B*C*D))

  # Values (equality constraints)
  bvec <- c(1, a, b, c, d, ab, bc, cd, da, ac, bd, abc, abd, acd, bcd)

  # Objective
  cvec <- with(X, A*B*C*D)

  # Eliminate unwanted constraints
  i <- which(!is.na(bvec))
  bvec <- bvec[i]
  X <- X[i]

  # Solve the max and min problems
  library(linprog)
  Amat <- t(as.matrix(X))
  dir <- rep("==", ncol(X)) # All constraints are equalities
  obj.min <- solveLP(cvec, bvec, Amat, const.dir = dir, lpSolve = TRUE, maximum = FALSE)
  obj.max <- solveLP(cvec, bvec, Amat, const.dir = dir, lpSolve = TRUE, maximum = TRUE)

  # Return the objects and the setup.
  list(Max = obj.max, Min = obj.min, X = X, bvec = bvec, cvec = cvec, Support = pm)
}
#
# An example.
#
Optima <- f(ab=1/2, bc=3/8, cd=1/4, da=1/8)
#
# Report the results.
#
for (s in c("Min", "Max")) {
  print(s)
  obj <- Optima[[s]]
  sol <- obj$solution
  if (is.null(sol)) {
print("*** Problem is infeasible. ***")
next
  }
  print(with(Optima$X, {
    c(`E[A]` = sum(A * sol),
      `E[B]` = sum(B * sol),
      `E[C]` = sum(C * sol),
      `E[D]` = sum(D * sol),
      `E[AB]` = sum(AB * sol),
      `E[BC]` = sum(BC * sol),
      `E[CD]` = sum(CD * sol),
      `E[DA]` = sum(DA * sol),
      `E[ABCD]` = sum(A*B*C*D * sol))
  }))
  #
  # Display the full distribution.
  #
  pm <- Optima$Support
  k <- length(pm)
  print(array(sol, dim = list(A=k, B=k, C=k, D=k),
         dimnames = list(A = pm, B = pm, C = pm, D = pm)))
}

