The heterogeneous autoregressive (HAR) model is a simple linear model that is commonly used to explain long memory in the realized volatility. However, as realized volatility has more complicated features such as conditional heteroscedasticity, leverage e ffect, and volatility clustering, it is necessary to extend the simple HAR model. Therefore, to better incorporate the stylized facts, we propose a threshold HAR model with GARCH errors, namely the THAR-GARCH model. That is, the THAR-GARCH model is a nonlinear model whose coe fficients vary according to a threshold value, and the conditional heteroscedasticity is explained through the GARCH errors. Model parameters are estimated using an iterative weighted least squares estimation method. Our simulation study supports the consistency of the iterative estimation method. In addition, we show that the proposed THAR-GARCH model has better forecasting power by applying to the realized volatility of major 21 stock indices around the world.