We use recent 36 observational Hubble data (OHD) in the redshift range .07 ≤ z ≤ 2.36, the joint light curves (JLA) sample, comprised of 740 type Ia supernovae (SNIa) in the redshift range 0.01 ≤ z ≤ 1.30, and their joint combination datasets to constrain anisotropic Bianchi type I (BI) dark energy model(DE). To estimate model parameters, we used Metropolis-Hasting algorithm to perform Monte Carlo Markov Chain analysis. We also compute the covariant matrix for BI dark energy model considering different datasets to compare the correlation between parameters of the model. To check the acceptability of our fittings, all results are compared with those obtained from 9 year WMAP as well as Planck (2015) collaboration. Our estimations show that at 68% confidence level, the dark energy equation of state (EOS) parameter for JLA data varies in quintessence-phantom region while for OHD and the joint combination of datasets only varies in phantom region. It is found that the current cosmic anisotropy is of order ∼ 10 −3 which imply that OHD and JLA datasets do not put tight constrain on this parameter. The deceleration parameter is obtained as q = −0.46 +0.89+0.36 −0.41−0.37 , q = −0.619 +0.12+0.20 −0.0.095−0.24 , and q = −0.52 +0.080+0.014 −0.046−0.15 for H(z), SNIa, and H(z)+SNIa data respectively.