-
Notifications
You must be signed in to change notification settings - Fork 269
Add Markov Switching GARCH volatility model #791
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Conversation
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #791 +/- ##
===========================================
- Coverage 99.70% 88.29% -11.41%
===========================================
Files 76 76
Lines 15740 15981 +241
Branches 1288 1303 +15
===========================================
- Hits 15694 14111 -1583
- Misses 17 1797 +1780
- Partials 29 73 +44
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for the good start. Should not use isinstance(..., MSGARCH)
and need to abstract to avoid this pattern.
args = (sigma2, backcast, var_bounds) | ||
ineq_constraints = constraint(a, b) | ||
if isinstance(self.volatility, MSGARCH): | ||
func = self.volatility.msgarch_loglikelihood # MS GARCH |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You should overrigt the default likelihood in the MSGARCH subclass so that you don't need to use a pattern like this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I understand the goal of avoiding isinstance checks for a cleaner design. The challenge is that the loglikelihood method resides in the ARCHModel class, while the new MSGARCH vol process is a subclass of VolatilityProcess in a separate module. This separation makes directly overriding the method difficult without significant architectural changes. To properly move the likelihood logic into the volatility classes would require refactoring much of the existing code, and will likely cause issues in other areas.
Sorry if im missing something obvious here, could you please provide some guidance on the best approach for this scenario?
|
||
from scipy.optimize import minimize | ||
else: | ||
func = self._loglikelihood # standard GARCH |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This isn't standard GARCH - it is every volatility process.
arch/univariate/base.py
Outdated
@@ -835,7 +845,10 @@ def fit( | |||
mp, vp, dp = self._parse_parameters(params) | |||
|
|||
resids = self.resids(mp) | |||
vol = np.zeros(resids.shape[0], dtype=float) | |||
if isinstance(self.volatility, MSGARCH): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This should be a new class-specific method. SOmething like
def _initialize_vol(self, t: int) -> Float64Array:
return np.zeros(t, dtype=float) # usual case or in the case of MSGARCH, a special shape
Obviously all use the same except MSGARCH (for now)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
implemented in "39ff059", added _initialise_vol method in VolatilityProcess with an overridden version in MSGARCH
arch/univariate/base.py
Outdated
vol_final = np.full(self._y.shape, np.nan, dtype=float) | ||
vol_final[first_obs:last_obs] = vol | ||
|
||
if isinstance(self.volatility, MSGARCH): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Need to abstract this somehow. Perhaps make it an attribute of a result. In the usual case, it would be a vector os 1s since there is only 1 state. Docstring would also explain in the usual case that this isn't meaningful unless the model has a MS structure.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In ‘a8f78d2’, I implemented a polymorphic compute_filtered_probs method in the VolatilityProcess base class:
• Default behaviour returns a vector of ones for single regime models.
• MSGARCH override implements the Hamilton filter, returning regime probabilities for each time step.
To calculate filtered probabilities for MSGARCH, I duplicated the logic from the MSGARCH loglikelihood function. This ensures that the filtered probabilities are computed using the same optimal parameters found during fitting, all without requiring isinstance checks in fit().
I duplicated the logic from the MSGARCH log likelihood function to ensure compute_filtered_probs works correctly first. Once this is confirmed, I plan to refactor by extracting the core Hamilton filter computation so it can be shared between the log likelihood and compute_filtered_probs.
I’ll be away for the next week, but I do plan on finishing this PR.
@@ -34,7 +34,7 @@ | |||
Literal, | |||
) | |||
from arch.univariate.distribution import Distribution, Normal | |||
from arch.univariate.volatility import ConstantVariance, VolatilityProcess | |||
from arch.univariate.volatility import ConstantVariance, VolatilityProcess, MSGARCH |
Check notice
Code scanning / CodeQL
Cyclic import
@@ -12,7 +12,8 @@ | |||
|
|||
import numpy as np | |||
from numpy.random import RandomState | |||
from scipy.special import gammaln | |||
from scipy.special import gammaln, logsumexp | |||
from scipy.stats import norm |
Check notice
Code scanning / CodeQL
Unused import
self.power = 2.0 # fixed for now | ||
|
||
self.num_params_single = 1 + self.p + self.q # parameters in a single regime | ||
self._num_params = self.num_params_single * self.k + self.k # regime specifc + transition matrix |
Check warning
Code scanning / CodeQL
Overwriting attribute in super-class or sub-class
# Initialise constant transition matrix | ||
self.transition_matrix = np.full((self.k, self.k), 1.0 / self.k) | ||
|
||
self._name = self._generate_name() |
Check warning
Code scanning / CodeQL
Overwriting attribute in super-class or sub-class
|
||
|
||
def msgarch_loglikelihood(self, parameters, resids, sigma2, backcast, var_bounds): | ||
from arch.univariate.base import _callback_info |
Check notice
Code scanning / CodeQL
Cyclic import
arch/univariate/volatility.py
Outdated
|
||
# GARCH params per regime | ||
garch_block = vp[:self.k*self.num_params_single] | ||
garch_params = garch_block.reshape(self.k, self.num_params_single) |
Check notice
Code scanning / CodeQL
Unused local variable
…CH. store filtered_probs in ARCHModelResult
14a9e23
to
a8f78d2
Compare
Summary:
Introduces MS GARCH to arch.univariate, supporting multiple regimes with separate GARCH(p,q) processes per regime.
Implementation:
The MS-GARCH model is implemented to closely match the structure and API of other Volatility Process models.
It supports 2 regimes and is estimated via Maximum Likelihood Estimation (MLE), without MCMC fitting (for now).
Starting values for the MLE are obtained via a gaussian mixture approximation of the returns. The loglikelihood is computed using the Hamilton filter to account for regime probabilities.
Conditional variances at each time step are calculated as the probability weighted average of the per regime GARCH variances.
One step ahead forecasting is implemented through the analytic forecast function.
All model parameters are constrained to ensure positivity and stationarity where apt.
Testing:
I added some basic tests to make sure the MS GARCH model runs, parameters look sensible, forecasts and variances are positive, and filtered probabilities sum to 1.
Definitely open to suggestions on extra checks or better ways to test it.
Possible Improvements:
During testing, the optimiser often raises a ConvergenceWarning (code 4: inequality constraints incompatible).Despite this warning, the fitted parameters and log likelihoods appear reasonable and don’t violate the constraints. I’ve investigated but haven’t pinpointed the cause, I’m open to suggestions on how to address or better understand this issue.
Implementing Bayesian estimation for more robust inference.
Currently limited to 2 regimes with p=q=1 chosen to provide a simple baseline model that integrates cleanly and is easy to test. This could be extended to more regimes, higher order GARCH processes, or non zero mean terms in the future.
Forecasting currently returns a weighted average of the conditional volatilities across regimes. Optionally, it could allow forecasts conditional on a single regime.
The GARCH recursion is the main computational bottleneck. performance could be improved with e.g. Numba or Cython.
This is my first contribution to the arch library so I’m happy to receive feedback. Thanks for maintaining such a great library, its been a great way to improve my knowledge of financial econometrics + coding! Let me know if there are any issues or you would need more for a first working version of MS GARCH functionality.