Estimation of covariate-dependent conditional covariance matrix in a high-dimensional space poses a challenge to contemporary statistical research. The existing kernel estimators may not be locally adaptive due to using a single bandwidth to explore the smoothness of all entries of the target matrix function. Moreover, the corresponding theory holds only for i.i.d. samples although in most of applications, the samples are dependent. In this paper, we propose a novel estimation scheme to overcome these obstacles by using techniques of factorization, thresholding and optimal shrinkage. Under certain regularity conditions, we show that the proposed estimator is consistent with the underlying matrix even when the sample is dependent. We conduct a set of simulation studies to show that the proposed estimator significantly outperforms its competitors. We apply the proposed procedure to the analysis of an asset return dataset, identifying a number of interesting volatility and co-volatility patterns across different time periods.