We investigate the problem of bounding counterfactual queries from an arbitrary collection of observational and experimental distributions and qualitative knowledge about the underlying data-generating model represented in the form of a causal diagram. We show that all counterfactual distributions in an arbitrary structural causal model (SCM) with finite discrete endogenous variables could be generated by a family of SCMs with the same causal diagram where unobserved (exogenous) variables are discrete with a finite domain. Utilizing this family of universal SCMs, we translate the problem of bounding counterfactuals into that of polynomial programming whose solution provides optimal bounds for the counterfactual query. Solving such polynomial programs is in general computationally expensive. We develop effective Monte Carlo algorithms to approximate optimal bounds from a combination of observational and experimental data.