<reponame>DanielBok/muarch
import os
from concurrent.futures import ProcessPoolExecutor
from functools import partial
from typing import Collection, Optional, Union

import numpy as np
import pandas as pd

from .exceptions import EndogInputError, InvalidModelError
from .summary import SummaryList
from .typings import Endog, Exog, RngGen
from .uarch import UArch


class MUArch:
    """
    Multi-univariate ARCH model. Unlike a multivariate ARCH model, this fits each univariate time series individually.
    Any simulations returns simulations of each univariate series column bound together.
    """

    def __init__(self, n: Union[int, Collection[UArch]], mean='Constant', lags=0, vol='GARCH', p=1, o=0, q=1, power=2.0,
                 dist='Normal', hold_back=None, scale=1):
        r"""
        Initializes the MUArch model.

        The MUArch model holds multiple univariate models which are determined during fitting. If the models are not
        specified, the global default options will be used. Models can be individually specified after initializing
        the MUArch instance.

        Parameters
        ----------
        n:  int or list of UArch models
            Number of univariate models to fit. Alternatively, a list of UArch (univariate) models can be specified.

        mean: { 'zero', 'constant', 'harx', 'har', 'ar', 'arx', 'ls' }, optional
            Name of the global default mean model.  Currently supported options are:

            * **Constant** (default) - Constant mean model
            * **Zero** - Zero mean model
            * **AR** - Autoregression model
            * **ARX** - Autoregression model with exogenous regressors. Falls back to **AR** if no exogenous regressors
            * **HAR** - Heterogeneous Autoregression model
            * **HARX** - Heterogeneous Autoregressions with exogenous regressors
            * **LS** - Least squares model

            For more information on the different models, check out the documentation at
            https://arch.readthedocs.io/en/latest/univariate/mean.html

        lags: int or list (int), optional
            Global default lag. Either a scalar integer value indicating lag length or a list of integers specifying
            lag locations.

        vol: { 'GARCH', 'ARCH', 'EGARCH', 'FIGARCH' and 'HARCH', 'CONSTANT' }, optional
            Name of the global default volatility model.  Currently supported options are:

            * **GARCH** (default) - Standard GARCH process which can be used to specify the following models:

                * ARCH(p)
                * GARCH(p,q)
                * GJR-GARCH(p,o,q)
                * AVARCH(p)
                * AVGARCH(p,q)
                * TARCH(p,o,q)
                * Models with arbitrary, pre-specified powers

            * **ARCH** - ARCH process
            * **EGARCH** - EGARCH process
            * **FIGARCH** - Fractionally Integrated (FI) GARCH process
            * **HARCH** - Heterogeneous ARCH process
            * **Constant** (default) - Constant volatility process

        p: int, optional
            Global default lag order of the symmetric innovation

        o: int, optional
            Global default lag order of the asymmetric innovation

        q: int, optional
            Global default lag order of lagged volatility or equivalent

        power: float, optional
            Global default power to use with the innovations, Default is 2.0, which produces ARCH and related
            models. Using 1.0 produces AVARCH and related models. Other powers can be specified, although these should
            be strictly positive, and usually larger than 0.25.

        dist:  { 'normal', 'gaussian', 'studentst', 't', 'skewstudent', 'skewt', 'ged', 'generalized error' }, optional
            Name of the global default distribution for the innovations.  Currently supported options are:

            * **normal**, **gaussian** (default) - Standard Normal distribution
            * **t**, **studentst** - Standardized Student’s distribution
            * **skewstudent**, **skewt** - Standardized Skewed Student’s distribution.
            * **ged**, **generalized error" - Generalized Error Distribution

        hold_back: int
            Global default. Number of observations at the start of the sample to exclude when estimating model
            parameters. Used when comparing models with different lag lengths to estimate on the common sample.

        scale: float
            Global default factor to scale data up or down by. This is useful when your data is too small leading to
            numerical errors when fitting. It will be used to scale simulation data
        """
        if isinstance(n, Collection):
            self._models = list(n)

            assert len(self) != 0, "If passing in a list of UArch models, list cannot be empty!"

            for i, m in enumerate(self._models):
                if not isinstance(m, UArch):
                    raise InvalidModelError(i)

            self._model_names = [str(i) for i in range(len(self))]

        elif isinstance(n, int):
            self._models = [UArch(mean, lags, vol, p, o, q, power, dist, hold_back, scale) for _ in range(n)]
            self._model_names = [str(i) for i in range(n)]
        else:
            raise TypeError('`n` must either be an integer specifying the number of models or a list of UArch models')

    def fit(self, y: Endog, x: Exog = None, update_freq=1, disp='off', cov_type='robust',
            show_warning=True, tol: float = None, options=None):
        """
        Fits the MUArch model.

        If finer control over the MUArch models is required, set the UArch models separately. Otherwise, method will
        set the default parameters.

        Parameters
        ----------
        y: {ndarray, Series}
            The dependent variable. If a vector is passed in, it is assumed that the same vector (endog) is used for
            all models. Otherwise, the last value of the shape must match the number of models

        x: list of {ndarray, None}, optional
            Exogenous regressors.  Ignored if model does not permit exogenous regressors. If passed in, the first shape
            must match the number of models.

        update_freq: int, optional
            Frequency of iteration updates.  Output is generated every `update_freq` iterations. Set to 0 to disable
            iterative output

        disp: 'final' or 'off' (default)
            Either 'final' to print optimization result or 'off' to display nothing

        cov_type: str, optional
            Estimation method of parameter covariance.  Supported options are 'robust', which does not assume the
            Information Matrix Equality holds and 'classic' which does. In the ARCH literature, 'robust' corresponds
            to Bollerslev-Wooldridge covariance estimator.

        show_warning: bool, optional
            Flag indicating whether convergence warnings should be shown.

        tol: float, optional
            Tolerance for termination

        options: dict, optional
            Options to pass to `scipy.optimize.minimize`.  Valid entries include 'ftol', 'eps', 'disp', and 'maxiter'

        Returns
        -------
        MUArch
            Fitted self instance
        """
        if y.ndim == 1:
            y = np.tile(y, len(self))

        if y.shape[1] != len(self) or y.ndim > 2:  # checks on input data
            raise EndogInputError

        if isinstance(y, pd.DataFrame):  # setting names for params later
            self._model_names = list(y.columns)
            y = y.values
        else:
            self._model_names = [str(i) for i in range(len(self))]

        for i in range(len(self)):
            yy = y[:, i]
            xx = x[i] if x is not None else None
            if xx is not None:
                xx = np.asarray(xx)
                if xx.ndim == 1:
                    xx = xx.reshape(len(xx), 1)

            self[i].fit(yy, xx, update_freq=update_freq, disp=disp, cov_type=cov_type, show_warning=show_warning,
                        tol=tol,
                        options=options)

    @property
    def params(self):
        return pd.DataFrame([m.params for m in self], self._model_names)

    def residuals(self, standardize=True) -> np.ndarray:
        """
        Model residuals

        The residuals will be burnt by the maximum lag of the underlying models. For example, given 3 models - AR(1),
        AR(10), Constant with 400 data points each, the residuals will be 399, 390 and 400 long. The function will
        cut off the first 10 data points in this instance.

        Parameters
        ----------
        standardize: bool, optional
            Whether to standardize residuals. Residuals are standardized by dividing it with the conditional volatility

        Returns
        -------
        ndarray
            Residuals
        """
        residuals = [mm.residuals(standardize) for mm in self]
        length = min([len(r) for r in residuals])

        return np.asarray([r[-length:] for r in residuals]).T

    def simulate(self,
                 nobs,
                 burn=500,
                 initial_value=None,
                 x=None,
                 initial_value_vol=None,
                 data_only=True,
                 custom_dist: Optional[Union[RngGen, np.ndarray]] = None):
        """
        Simulates data from the multiple ARMA-GARCH models

        Parameters
        ----------
        nobs: int
            Length of series to simulate

        burn: int, optional
            Number of values to simulate to initialize the model and remove dependence on initial values

        initial_value: {ndarray, float}, optional
            Either a scalar value or `max(lags)` array set of initial values to use when initializing the model.
            If omitted, 0.0 is used. If array, the last column must be of the same size as the number of models

        x: {ndarray, list of ndarray}, optional
            If supplied as a list, this list should have the same number of elements as the number of models in the
            MUArch model. Each array inside is the specified exogenous variable for that  particular model and
            this must be a nobs + burn by k matrix of exogenous variables to include in the simulation. Otherwise,
            leave the value as `None` to indicate no exogenous variables are used for simulation in the model.

            If an array is supplied directly, it means every model has an exogenous variable associated with it. In
            this case, it should be a 3 dimensional tensor where the first dimension represents the number of models.

        initial_value_vol: {ndarray, float}, optional
            An array or scalar to use when initializing the volatility process. If array, the last column must be of the
            same size as the number of models

        data_only: bool, default True
            If True, this returns only the simulated data, omits the volatility and error. In this case, it will return
            as a numpy array. Otherwise, it returns a data frame with the data, volatility and error

        custom_dist: {ndarray, Callable}, optional
            Optional density from which to simulate the innovations (Distribution) in the GARCH models. This is useful
            when working with the copula-GARCH model where each univariate model innovations has dependence on others.
            It is assumed that the values supplied are standardized [0, 1] innovations instead of the unstandardized
            residuals.

            The shape of the array must be at least as long as the simulation size required after accounting for burn
            and type of innovation process. If unsure, use :code:`simulation_size_required` to check. It must also
            account for the number of dimensions of the MUArch model. For example, if MUArch model is simulating a
            horizon of 120 time steps, 10000 trials and has 5 UArch models, the shape of the numpy array should be
            (120, 10000, 5).

            If a random number generator function is passed in, ensure that it only takes only argument and returns
            a numpy array. The argument can be an integer or a tuple of integers. In this case, the size will be
            automatically derived to save the user the trouble.

        Returns
        -------
        simulated_data : {List[DataFrame], ndarray}
            List of DataFrame with columns data containing the simulated values, volatility, containing the conditional
            volatility and errors containing the errors used in the simulation

        See Also
        --------
        UArch.simulation_horizon_required: Calculates the simulation size required
        """

        initial_value, initial_value_vol, x, custom_dist = self._format_simulation_parameters(nobs,
                                                                                              burn,
                                                                                              None,
                                                                                              initial_value,
                                                                                              initial_value_vol,
                                                                                              x,
                                                                                              custom_dist)

        sims = []
        for i in range(len(self)):
            r = custom_dist[..., i] if custom_dist is not None else None
            xx = x[i]
            iv = initial_value[i]
            ivv = initial_value_vol[i]

            sims.append(self[i].simulate(nobs, burn, iv, xx, ivv, data_only, None, r))

        if data_only:
            sims = np.asarray(sims).T

        return sims

    def simulate_mc(self,
                    nobs,
                    reps,
                    burn=500,
                    initial_value=None,
                    x=None,
                    initial_value_vol=None,
                    custom_dist: Optional[Union[RngGen, np.ndarray]] = None,
                    n_jobs: Optional[int] = None):
        """
        Simulates data from the multiple ARCH-GARCH models.

        This function is specially crafted for Monte-Carlo simulations.

        Parameters
        ----------
        nobs: int
            Length of series to simulate

        reps: int
            Number of repetitions

        burn: int, optional
            Number of values to simulate to initialize the model and remove dependence on initial values

        initial_value: {ndarray, float}, optional
            Either a scalar value or `max(lags)` array set of initial values to use when initializing the model.
            If omitted, 0.0 is used. If array, the last column must be of the same size as the number of models

        x: {ndarray, list of ndarray}, optional
            If supplied as a list, this list should have the same number of elements as the number of models in the
            MUArch model. Each array inside is the specified exogenous variable for that  particular model and
            this must be a nobs + burn by k matrix of exogenous variables to include in the simulation. Otherwise,
            leave the value as `None` to indicate no exogenous variables are used for simulation in the model.

            If an array is supplied directly, it means every model has an exogenous variable associated with it. In
            this case, it should be a 3 dimensional tensor where the first dimension represents the number of models.

        initial_value_vol: {ndarray, float}, optional
            An array or scalar to use when initializing the volatility process. If array, the last column must be of the
            same size as the number of models

        custom_dist: {ndarray, Callable}, optional
            Optional density from which to simulate the innovations (Distribution) in the GARCH models. This is useful
            when working with the copula-GARCH model where each univariate model innovations has dependence on others.
            It is assumed that the values supplied are standardized [0, 1] innovations instead of the unstandardized
            residuals.

            The shape of the array must be at least as long as the simulation size required after accounting for burn
            and type of innovation process. If unsure, use :code:`simulation_size_required` to check.

            If a random number generator function is passed in, he size will be automatically derived to save the user
            the trouble. However, the function must:

            * take as it first argument an integer or a tuple of integer
            * have other parameters that are optional
            * return a numpy array

        n_jobs: int or None, optional
            The number of jobs to run in parallel for simulation. This is particularly useful when simulating large
            number of repetitions with more than 1 dimension. None defaults to using 1 processor. Any numbers less or
            equal to 0 means to use all processors. Even if a large number is used, it will be capped at the
            maximum number of processors available.


        Returns
        -------
        simulated_data : numpy array
            Array containing simulated data from the Monte Carlo Simulation
        """
        if not isinstance(reps, int) or reps < 1:
            raise ValueError('reps must be an integer greater than 0')

        # setting up job processes
        if n_jobs is None:
            n_jobs = 1
        elif not isinstance(n_jobs, int):
            raise ValueError('`n_jobs` must be an integer or None')
        else:
            n_jobs = min(os.cpu_count(), os.cpu_count() if n_jobs <= 0 else n_jobs)

        initial_value, initial_value_vol, x, custom_dist = self._format_simulation_parameters(nobs,
                                                                                              burn,
                                                                                              reps,
                                                                                              initial_value,
                                                                                              initial_value_vol,
                                                                                              x,
                                                                                              custom_dist)

        # formulate function calls to get simulated data
        functions = []
        for i in range(len(self)):
            r = custom_dist[..., i] if custom_dist is not None else None
            xx = x[i]
            iv = initial_value[i]
            ivv = initial_value_vol[i]

            func = partial(self[i].simulate_mc, nobs, reps, burn, iv, xx, ivv, None, r)
            functions.append(func)

        # switch between sequential and multiprocessing
        if n_jobs == 1:  # sequential
            sims = np.zeros((nobs, reps, len(self)))
            for i, func in enumerate(functions):
                sims[:, :, i] = func()

            return sims
        else:  # multiprocessing
            results = []

            with ProcessPoolExecutor(n_jobs) as P:
                for func in functions:
                    results.append(P.submit(func))

            results = np.asarray([r.result() for r in results])
            return np.moveaxis(results, 0, -1)

    def summary(self, short=False, dp=4):
        """
        Summary of fitted models

        :param short: bool, default False
            Whether to show short summary or full summary.
        :param dp: int, default 4
            Number of decimal places to show in short summary
        :return: SummaryList
            summary of fitted models
        """
        if short:
            return pd.DataFrame([m.summary(short, dp) for m in self], self._model_names)
        else:
            smry = SummaryList()
            for i, name in enumerate(self._model_names):
                s = self[i].summary(short, dp)
                s.add_header(name)
                smry.append(s)
            return smry

    def _format_simulation_parameters(self,
                                      nobs,
                                      burn,
                                      reps,
                                      initial_value,
                                      initial_value_vol,
                                      x,
                                      custom_dist):
        """Ensures the simulation are okay and return the formatted parameters"""

        # Format initial value
        if initial_value is None or isinstance(initial_value, (float, int)):
            initial_value = np.asarray([initial_value] * len(self))

        # Format initial volatility value
        if initial_value_vol is None or isinstance(initial_value, (float, int)):
            initial_value_vol = np.asarray([initial_value_vol] * len(self))

        # Format exogenous regressors
        if x is not None:
            if isinstance(x, np.ndarray):
                assert x.ndim == 3, ("If numpy array passed in as MUArch `x` (exog), make sure it is a 3 "
                                     "dimensional array where the last dimension is the number of elements")
                x = [x[i] for i in range(len(self))]
            else:
                x = list(x)

                assert len(x) == len(self), ("Exogenous variable's list should have the same number of elements as the "
                                             "number of models")

                for i, xx in enumerate(x):
                    if xx is None:
                        continue
                    elif isinstance(x[0], pd.DataFrame):
                        # convert all data frames to numpy array then move the model axis (first) to last
                        xx = np.moveaxis(np.array([df.values for df in x]), 0, -1)
                    else:
                        xx = np.asarray(xx)

                    if xx.ndim == 1:
                        xx = xx.reshape(-1, 1)
                    x[i] = xx
        else:
            x = [None] * len(self)

        # Format custom distribution
        if custom_dist is not None:
            horizon = max(m.simulation_horizon(nobs, burn) for m in self)
            if isinstance(reps, int):  # for monte carlo
                size = horizon, reps
            else:
                size = horizon

            if callable(custom_dist):
                custom_dist = np.asarray(custom_dist(size))
            else:
                custom_dist = np.asarray(custom_dist)

            if custom_dist.ndim == 1:
                custom_dist = custom_dist.reshape(1, -1)

            dim = custom_dist.shape[-1]
            assert dim == len(self), (f"expected generator to return array with {len(self)} vectors in the last axis "
                                      f"(of shape) but got {dim}")

        return initial_value, initial_value_vol, x, custom_dist

    def _repr_html_(self):
        header = '<h2>Arch Models</h2>'
        models = '<br/>'.join(m._repr_html_() for m in self)
        return header + models

    def __getitem__(self, i: int):
        return self._models[i]

    def __setitem__(self, i: int, model: UArch):
        assert isinstance(model, UArch), 'only UArch models are allowed to be set as a component of the MUArch object'
        self._models[i] = model

    def __repr__(self):
        return self.__str__() + f', {hex(id(self))}'

    def __str__(self):
        txt = f'{"Arch Models":^80}\n' + '=' * 80 + '\n\n'
        txt += '\n\n\n'.join(str(m) for m in self)
        return txt

    def __len__(self):
        return len(self._models)

import os
from concurrent.futures import ProcessPoolExecutor
from functools import partial
from typing import Collection, Optional, Union

import numpy as np
import pandas as pd

from .exceptions import EndogInputError, InvalidModelError
from .summary import SummaryList
from .typings import Endog, Exog, RngGen
from .uarch import UArch


class MUArch:
    """
    Multi-univariate ARCH model. Unlike a multivariate ARCH model, this fits each univariate time series individually.
    Any simulations returns simulations of each univariate series column bound together.
    """

    def __init__(self, n: Union[int, Collection[UArch]], mean='Constant', lags=0, vol='GARCH', p=1, o=0, q=1, power=2.0,
                 dist='Normal', hold_back=None, scale=1):
        r"""
        Initializes the MUArch model.

        The MUArch model holds multiple univariate models which are determined during fitting. If the models are not
        specified, the global default options will be used. Models can be individually specified after initializing
        the MUArch instance.

        Parameters
        ----------
        n:  int or list of UArch models
            Number of univariate models to fit. Alternatively, a list of UArch (univariate) models can be specified.

        mean: { 'zero', 'constant', 'harx', 'har', 'ar', 'arx', 'ls' }, optional
            Name of the global default mean model.  Currently supported options are:

            * **Constant** (default) - Constant mean model
            * **Zero** - Zero mean model
            * **AR** - Autoregression model
            * **ARX** - Autoregression model with exogenous regressors. Falls back to **AR** if no exogenous regressors
            * **HAR** - Heterogeneous Autoregression model
            * **HARX** - Heterogeneous Autoregressions with exogenous regressors
            * **LS** - Least squares model

            For more information on the different models, check out the documentation at
            https://arch.readthedocs.io/en/latest/univariate/mean.html

        lags: int or list (int), optional
            Global default lag. Either a scalar integer value indicating lag length or a list of integers specifying
            lag locations.

        vol: { 'GARCH', 'ARCH', 'EGARCH', 'FIGARCH' and 'HARCH', 'CONSTANT' }, optional
            Name of the global default volatility model.  Currently supported options are:

            * **GARCH** (default) - Standard GARCH process which can be used to specify the following models:

                * ARCH(p)
                * GARCH(p,q)
                * GJR-GARCH(p,o,q)
                * AVARCH(p)
                * AVGARCH(p,q)
                * TARCH(p,o,q)
                * Models with arbitrary, pre-specified powers

            * **ARCH** - ARCH process
            * **EGARCH** - EGARCH process
            * **FIGARCH** - Fractionally Integrated (FI) GARCH process
            * **HARCH** - Heterogeneous ARCH process
            * **Constant** (default) - Constant volatility process

        p: int, optional
            Global default lag order of the symmetric innovation

        o: int, optional
            Global default lag order of the asymmetric innovation

        q: int, optional
            Global default lag order of lagged volatility or equivalent

        power: float, optional
            Global default power to use with the innovations, Default is 2.0, which produces ARCH and related
            models. Using 1.0 produces AVARCH and related models. Other powers can be specified, although these should
            be strictly positive, and usually larger than 0.25.

        dist:  { 'normal', 'gaussian', 'studentst', 't', 'skewstudent', 'skewt', 'ged', 'generalized error' }, optional
            Name of the global default distribution for the innovations.  Currently supported options are:

            * **normal**, **gaussian** (default) - Standard Normal distribution
            * **t**, **studentst** - Standardized Student’s distribution
            * **skewstudent**, **skewt** - Standardized Skewed Student’s distribution.
            * **ged**, **generalized error" - Generalized Error Distribution

        hold_back: int
            Global default. Number of observations at the start of the sample to exclude when estimating model
            parameters. Used when comparing models with different lag lengths to estimate on the common sample.

        scale: float
            Global default factor to scale data up or down by. This is useful when your data is too small leading to
            numerical errors when fitting. It will be used to scale simulation data
        """
        if isinstance(n, Collection):
            self._models = list(n)

            assert len(self) != 0, "If passing in a list of UArch models, list cannot be empty!"

            for i, m in enumerate(self._models):
                if not isinstance(m, UArch):
                    raise InvalidModelError(i)

            self._model_names = [str(i) for i in range(len(self))]

        elif isinstance(n, int):
            self._models = [UArch(mean, lags, vol, p, o, q, power, dist, hold_back, scale) for _ in range(n)]
            self._model_names = [str(i) for i in range(n)]
        else:
            raise TypeError('`n` must either be an integer specifying the number of models or a list of UArch models')

    def fit(self, y: Endog, x: Exog = None, update_freq=1, disp='off', cov_type='robust',
            show_warning=True, tol: float = None, options=None):
        """
        Fits the MUArch model.

        If finer control over the MUArch models is required, set the UArch models separately. Otherwise, method will
        set the default parameters.

        Parameters
        ----------
        y: {ndarray, Series}
            The dependent variable. If a vector is passed in, it is assumed that the same vector (endog) is used for
            all models. Otherwise, the last value of the shape must match the number of models

        x: list of {ndarray, None}, optional
            Exogenous regressors.  Ignored if model does not permit exogenous regressors. If passed in, the first shape
            must match the number of models.

        update_freq: int, optional
            Frequency of iteration updates.  Output is generated every `update_freq` iterations. Set to 0 to disable
            iterative output

        disp: 'final' or 'off' (default)
            Either 'final' to print optimization result or 'off' to display nothing

        cov_type: str, optional
            Estimation method of parameter covariance.  Supported options are 'robust', which does not assume the
            Information Matrix Equality holds and 'classic' which does. In the ARCH literature, 'robust' corresponds
            to Bollerslev-Wooldridge covariance estimator.

        show_warning: bool, optional
            Flag indicating whether convergence warnings should be shown.

        tol: float, optional
            Tolerance for termination

        options: dict, optional
            Options to pass to `scipy.optimize.minimize`.  Valid entries include 'ftol', 'eps', 'disp', and 'maxiter'

        Returns
        -------
        MUArch
            Fitted self instance
        """
        if y.ndim == 1:
            y = np.tile(y, len(self))

        if y.shape[1] != len(self) or y.ndim > 2:  # checks on input data
            raise EndogInputError

        if isinstance(y, pd.DataFrame):  # setting names for params later
            self._model_names = list(y.columns)
            y = y.values
        else:
            self._model_names = [str(i) for i in range(len(self))]

        for i in range(len(self)):
            yy = y[:, i]
            xx = x[i] if x is not None else None
            if xx is not None:
                xx = np.asarray(xx)
                if xx.ndim == 1:
                    xx = xx.reshape(len(xx), 1)

            self[i].fit(yy, xx, update_freq=update_freq, disp=disp, cov_type=cov_type, show_warning=show_warning,
                        tol=tol,
                        options=options)

    @property
    def params(self):
        return pd.DataFrame([m.params for m in self], self._model_names)

    def residuals(self, standardize=True) -> np.ndarray:
        """
        Model residuals

        The residuals will be burnt by the maximum lag of the underlying models. For example, given 3 models - AR(1),
        AR(10), Constant with 400 data points each, the residuals will be 399, 390 and 400 long. The function will
        cut off the first 10 data points in this instance.

        Parameters
        ----------
        standardize: bool, optional
            Whether to standardize residuals. Residuals are standardized by dividing it with the conditional volatility

        Returns
        -------
        ndarray
            Residuals
        """
        residuals = [mm.residuals(standardize) for mm in self]
        length = min([len(r) for r in residuals])

        return np.asarray([r[-length:] for r in residuals]).T

    def simulate(self,
                 nobs,
                 burn=500,
                 initial_value=None,
                 x=None,
                 initial_value_vol=None,
                 data_only=True,
                 custom_dist: Optional[Union[RngGen, np.ndarray]] = None):
        """
        Simulates data from the multiple ARMA-GARCH models

        Parameters
        ----------
        nobs: int
            Length of series to simulate

        burn: int, optional
            Number of values to simulate to initialize the model and remove dependence on initial values

        initial_value: {ndarray, float}, optional
            Either a scalar value or `max(lags)` array set of initial values to use when initializing the model.
            If omitted, 0.0 is used. If array, the last column must be of the same size as the number of models

        x: {ndarray, list of ndarray}, optional
            If supplied as a list, this list should have the same number of elements as the number of models in the
            MUArch model. Each array inside is the specified exogenous variable for that  particular model and
            this must be a nobs + burn by k matrix of exogenous variables to include in the simulation. Otherwise,
            leave the value as `None` to indicate no exogenous variables are used for simulation in the model.

            If an array is supplied directly, it means every model has an exogenous variable associated with it. In
            this case, it should be a 3 dimensional tensor where the first dimension represents the number of models.

        initial_value_vol: {ndarray, float}, optional
            An array or scalar to use when initializing the volatility process. If array, the last column must be of the
            same size as the number of models

        data_only: bool, default True
            If True, this returns only the simulated data, omits the volatility and error. In this case, it will return
            as a numpy array. Otherwise, it returns a data frame with the data, volatility and error

        custom_dist: {ndarray, Callable}, optional
            Optional density from which to simulate the innovations (Distribution) in the GARCH models. This is useful
            when working with the copula-GARCH model where each univariate model innovations has dependence on others.
            It is assumed that the values supplied are standardized [0, 1] innovations instead of the unstandardized
            residuals.

            The shape of the array must be at least as long as the simulation size required after accounting for burn
            and type of innovation process. If unsure, use :code:`simulation_size_required` to check. It must also
            account for the number of dimensions of the MUArch model. For example, if MUArch model is simulating a
            horizon of 120 time steps, 10000 trials and has 5 UArch models, the shape of the numpy array should be
            (120, 10000, 5).

            If a random number generator function is passed in, ensure that it only takes only argument and returns
            a numpy array. The argument can be an integer or a tuple of integers. In this case, the size will be
            automatically derived to save the user the trouble.

        Returns
        -------
        simulated_data : {List[DataFrame], ndarray}
            List of DataFrame with columns data containing the simulated values, volatility, containing the conditional
            volatility and errors containing the errors used in the simulation

        See Also
        --------
        UArch.simulation_horizon_required: Calculates the simulation size required
        """

        initial_value, initial_value_vol, x, custom_dist = self._format_simulation_parameters(nobs,
                                                                                              burn,
                                                                                              None,
                                                                                              initial_value,
                                                                                              initial_value_vol,
                                                                                              x,
                                                                                              custom_dist)

        sims = []
        for i in range(len(self)):
            r = custom_dist[..., i] if custom_dist is not None else None
            xx = x[i]
            iv = initial_value[i]
            ivv = initial_value_vol[i]

            sims.append(self[i].simulate(nobs, burn, iv, xx, ivv, data_only, None, r))

        if data_only:
            sims = np.asarray(sims).T

        return sims

    def simulate_mc(self,
                    nobs,
                    reps,
                    burn=500,
                    initial_value=None,
                    x=None,
                    initial_value_vol=None,
                    custom_dist: Optional[Union[RngGen, np.ndarray]] = None,
                    n_jobs: Optional[int] = None):
        """
        Simulates data from the multiple ARCH-GARCH models.

        This function is specially crafted for Monte-Carlo simulations.

        Parameters
        ----------
        nobs: int
            Length of series to simulate

        reps: int
            Number of repetitions

        burn: int, optional
            Number of values to simulate to initialize the model and remove dependence on initial values

        initial_value: {ndarray, float}, optional
            Either a scalar value or `max(lags)` array set of initial values to use when initializing the model.
            If omitted, 0.0 is used. If array, the last column must be of the same size as the number of models

        x: {ndarray, list of ndarray}, optional
            If supplied as a list, this list should have the same number of elements as the number of models in the
            MUArch model. Each array inside is the specified exogenous variable for that  particular model and
            this must be a nobs + burn by k matrix of exogenous variables to include in the simulation. Otherwise,
            leave the value as `None` to indicate no exogenous variables are used for simulation in the model.

            If an array is supplied directly, it means every model has an exogenous variable associated with it. In
            this case, it should be a 3 dimensional tensor where the first dimension represents the number of models.

        initial_value_vol: {ndarray, float}, optional
            An array or scalar to use when initializing the volatility process. If array, the last column must be of the
            same size as the number of models

        custom_dist: {ndarray, Callable}, optional
            Optional density from which to simulate the innovations (Distribution) in the GARCH models. This is useful
            when working with the copula-GARCH model where each univariate model innovations has dependence on others.
            It is assumed that the values supplied are standardized [0, 1] innovations instead of the unstandardized
            residuals.

            The shape of the array must be at least as long as the simulation size required after accounting for burn
            and type of innovation process. If unsure, use :code:`simulation_size_required` to check.

            If a random number generator function is passed in, he size will be automatically derived to save the user
            the trouble. However, the function must:

            * take as it first argument an integer or a tuple of integer
            * have other parameters that are optional
            * return a numpy array

        n_jobs: int or None, optional
            The number of jobs to run in parallel for simulation. This is particularly useful when simulating large
            number of repetitions with more than 1 dimension. None defaults to using 1 processor. Any numbers less or
            equal to 0 means to use all processors. Even if a large number is used, it will be capped at the
            maximum number of processors available.


        Returns
        -------
        simulated_data : numpy array
            Array containing simulated data from the Monte Carlo Simulation
        """
        if not isinstance(reps, int) or reps < 1:
            raise ValueError('reps must be an integer greater than 0')

        # setting up job processes
        if n_jobs is None:
            n_jobs = 1
        elif not isinstance(n_jobs, int):
            raise ValueError('`n_jobs` must be an integer or None')
        else:
            n_jobs = min(os.cpu_count(), os.cpu_count() if n_jobs <= 0 else n_jobs)

        initial_value, initial_value_vol, x, custom_dist = self._format_simulation_parameters(nobs,
                                                                                              burn,
                                                                                              reps,
                                                                                              initial_value,
                                                                                              initial_value_vol,
                                                                                              x,
                                                                                              custom_dist)

        # formulate function calls to get simulated data
        functions = []
        for i in range(len(self)):
            r = custom_dist[..., i] if custom_dist is not None else None
            xx = x[i]
            iv = initial_value[i]
            ivv = initial_value_vol[i]

            func = partial(self[i].simulate_mc, nobs, reps, burn, iv, xx, ivv, None, r)
            functions.append(func)

        # switch between sequential and multiprocessing
        if n_jobs == 1:  # sequential
            sims = np.zeros((nobs, reps, len(self)))
            for i, func in enumerate(functions):
                sims[:, :, i] = func()

            return sims
        else:  # multiprocessing
            results = []

            with ProcessPoolExecutor(n_jobs) as P:
                for func in functions:
                    results.append(P.submit(func))

            results = np.asarray([r.result() for r in results])
            return np.moveaxis(results, 0, -1)

    def summary(self, short=False, dp=4):
        """
        Summary of fitted models

        :param short: bool, default False
            Whether to show short summary or full summary.
        :param dp: int, default 4
            Number of decimal places to show in short summary
        :return: SummaryList
            summary of fitted models
        """
        if short:
            return pd.DataFrame([m.summary(short, dp) for m in self], self._model_names)
        else:
            smry = SummaryList()
            for i, name in enumerate(self._model_names):
                s = self[i].summary(short, dp)
                s.add_header(name)
                smry.append(s)
            return smry

    def _format_simulation_parameters(self,
                                      nobs,
                                      burn,
                                      reps,
                                      initial_value,
                                      initial_value_vol,
                                      x,
                                      custom_dist):
        """Ensures the simulation are okay and return the formatted parameters"""

        # Format initial value
        if initial_value is None or isinstance(initial_value, (float, int)):
            initial_value = np.asarray([initial_value] * len(self))

        # Format initial volatility value
        if initial_value_vol is None or isinstance(initial_value, (float, int)):
            initial_value_vol = np.asarray([initial_value_vol] * len(self))

        # Format exogenous regressors
        if x is not None:
            if isinstance(x, np.ndarray):
                assert x.ndim == 3, ("If numpy array passed in as MUArch `x` (exog), make sure it is a 3 "
                                     "dimensional array where the last dimension is the number of elements")
                x = [x[i] for i in range(len(self))]
            else:
                x = list(x)

                assert len(x) == len(self), ("Exogenous variable's list should have the same number of elements as the "
                                             "number of models")

                for i, xx in enumerate(x):
                    if xx is None:
                        continue
                    elif isinstance(x[0], pd.DataFrame):
                        # convert all data frames to numpy array then move the model axis (first) to last
                        xx = np.moveaxis(np.array([df.values for df in x]), 0, -1)
                    else:
                        xx = np.asarray(xx)

                    if xx.ndim == 1:
                        xx = xx.reshape(-1, 1)
                    x[i] = xx
        else:
            x = [None] * len(self)

        # Format custom distribution
        if custom_dist is not None:
            horizon = max(m.simulation_horizon(nobs, burn) for m in self)
            if isinstance(reps, int):  # for monte carlo
                size = horizon, reps
            else:
                size = horizon

            if callable(custom_dist):
                custom_dist = np.asarray(custom_dist(size))
            else:
                custom_dist = np.asarray(custom_dist)

            if custom_dist.ndim == 1:
                custom_dist = custom_dist.reshape(1, -1)

            dim = custom_dist.shape[-1]
            assert dim == len(self), (f"expected generator to return array with {len(self)} vectors in the last axis "
                                      f"(of shape) but got {dim}")

        return initial_value, initial_value_vol, x, custom_dist

    def _repr_html_(self):
        header = '<h2>Arch Models</h2>'
        models = '<br/>'.join(m._repr_html_() for m in self)
        return header + models

    def __getitem__(self, i: int):
        return self._models[i]

    def __setitem__(self, i: int, model: UArch):
        assert isinstance(model, UArch), 'only UArch models are allowed to be set as a component of the MUArch object'
        self._models[i] = model

    def __repr__(self):
        return self.__str__() + f', {hex(id(self))}'

    def __str__(self):
        txt = f'{"Arch Models":^80}\n' + '=' * 80 + '\n\n'
        txt += '\n\n\n'.join(str(m) for m in self)
        return txt

    def __len__(self):
        return len(self._models)
Multi-univariate ARCH model. Unlike a multivariate ARCH model, this fits each univariate time series individually.     Any simulations returns simulations of each univariate series column bound together. Initializes the MUArch model.          The MUArch model holds multiple univariate models which are determined during fitting. If the models are not         specified, the global default options will be used. Models can be individually specified after initializing         the MUArch instance.          Parameters         ----------         n:  int or list of UArch models             Number of univariate models to fit. Alternatively, a list of UArch (univariate) models can be specified.          mean: { 'zero', 'constant', 'harx', 'har', 'ar', 'arx', 'ls' }, optional             Name of the global default mean model.  Currently supported options are:              * **Constant** (default) - Constant mean model             * **Zero** - Zero mean model             * **AR** - Autoregression model             * **ARX** - Autoregression model with exogenous regressors. Falls back to **AR** if no exogenous regressors             * **HAR** - Heterogeneous Autoregression model             * **HARX** - Heterogeneous Autoregressions with exogenous regressors             * **LS** - Least squares model              For more information on the different models, check out the documentation at             https://arch.readthedocs.io/en/latest/univariate/mean.html          lags: int or list (int), optional             Global default lag. Either a scalar integer value indicating lag length or a list of integers specifying             lag locations.          vol: { 'GARCH', 'ARCH', 'EGARCH', 'FIGARCH' and 'HARCH', 'CONSTANT' }, optional             Name of the global default volatility model.  Currently supported options are:              * **GARCH** (default) - Standard GARCH process which can be used to specify the following models:                  * ARCH(p)                 * GARCH(p,q)                 * GJR-GARCH(p,o,q)                 * AVARCH(p)                 * AVGARCH(p,q)                 * TARCH(p,o,q)                 * Models with arbitrary, pre-specified powers              * **ARCH** - ARCH process             * **EGARCH** - EGARCH process             * **FIGARCH** - Fractionally Integrated (FI) GARCH process             * **HARCH** - Heterogeneous ARCH process             * **Constant** (default) - Constant volatility process          p: int, optional             Global default lag order of the symmetric innovation          o: int, optional             Global default lag order of the asymmetric innovation          q: int, optional             Global default lag order of lagged volatility or equivalent          power: float, optional             Global default power to use with the innovations, Default is 2.0, which produces ARCH and related             models. Using 1.0 produces AVARCH and related models. Other powers can be specified, although these should             be strictly positive, and usually larger than 0.25.          dist:  { 'normal', 'gaussian', 'studentst', 't', 'skewstudent', 'skewt', 'ged', 'generalized error' }, optional             Name of the global default distribution for the innovations.  Currently supported options are:              * **normal**, **gaussian** (default) - Standard Normal distribution             * **t**, **studentst** - Standardized Student’s distribution             * **skewstudent**, **skewt** - Standardized Skewed Student’s distribution.             * **ged**, **generalized error" - Generalized Error Distribution          hold_back: int             Global default. Number of observations at the start of the sample to exclude when estimating model             parameters. Used when comparing models with different lag lengths to estimate on the common sample.          scale: float             Global default factor to scale data up or down by. This is useful when your data is too small leading to             numerical errors when fitting. It will be used to scale simulation data Fits the MUArch model.          If finer control over the MUArch models is required, set the UArch models separately. Otherwise, method will         set the default parameters.          Parameters         ----------         y: {ndarray, Series}             The dependent variable. If a vector is passed in, it is assumed that the same vector (endog) is used for             all models. Otherwise, the last value of the shape must match the number of models          x: list of {ndarray, None}, optional             Exogenous regressors.  Ignored if model does not permit exogenous regressors. If passed in, the first shape             must match the number of models.          update_freq: int, optional             Frequency of iteration updates.  Output is generated every `update_freq` iterations. Set to 0 to disable             iterative output          disp: 'final' or 'off' (default)             Either 'final' to print optimization result or 'off' to display nothing          cov_type: str, optional             Estimation method of parameter covariance.  Supported options are 'robust', which does not assume the             Information Matrix Equality holds and 'classic' which does. In the ARCH literature, 'robust' corresponds             to Bollerslev-Wooldridge covariance estimator.          show_warning: bool, optional             Flag indicating whether convergence warnings should be shown.          tol: float, optional             Tolerance for termination          options: dict, optional             Options to pass to `scipy.optimize.minimize`.  Valid entries include 'ftol', 'eps', 'disp', and 'maxiter'          Returns         -------         MUArch             Fitted self instance # checks on input data # setting names for params later Model residuals          The residuals will be burnt by the maximum lag of the underlying models. For example, given 3 models - AR(1),         AR(10), Constant with 400 data points each, the residuals will be 399, 390 and 400 long. The function will         cut off the first 10 data points in this instance.          Parameters         ----------         standardize: bool, optional             Whether to standardize residuals. Residuals are standardized by dividing it with the conditional volatility          Returns         -------         ndarray             Residuals Simulates data from the multiple ARMA-GARCH models          Parameters         ----------         nobs: int             Length of series to simulate          burn: int, optional             Number of values to simulate to initialize the model and remove dependence on initial values          initial_value: {ndarray, float}, optional             Either a scalar value or `max(lags)` array set of initial values to use when initializing the model.             If omitted, 0.0 is used. If array, the last column must be of the same size as the number of models          x: {ndarray, list of ndarray}, optional             If supplied as a list, this list should have the same number of elements as the number of models in the             MUArch model. Each array inside is the specified exogenous variable for that  particular model and             this must be a nobs + burn by k matrix of exogenous variables to include in the simulation. Otherwise,             leave the value as `None` to indicate no exogenous variables are used for simulation in the model.              If an array is supplied directly, it means every model has an exogenous variable associated with it. In             this case, it should be a 3 dimensional tensor where the first dimension represents the number of models.          initial_value_vol: {ndarray, float}, optional             An array or scalar to use when initializing the volatility process. If array, the last column must be of the             same size as the number of models          data_only: bool, default True             If True, this returns only the simulated data, omits the volatility and error. In this case, it will return             as a numpy array. Otherwise, it returns a data frame with the data, volatility and error          custom_dist: {ndarray, Callable}, optional             Optional density from which to simulate the innovations (Distribution) in the GARCH models. This is useful             when working with the copula-GARCH model where each univariate model innovations has dependence on others.             It is assumed that the values supplied are standardized [0, 1] innovations instead of the unstandardized             residuals.              The shape of the array must be at least as long as the simulation size required after accounting for burn             and type of innovation process. If unsure, use :code:`simulation_size_required` to check. It must also             account for the number of dimensions of the MUArch model. For example, if MUArch model is simulating a             horizon of 120 time steps, 10000 trials and has 5 UArch models, the shape of the numpy array should be             (120, 10000, 5).              If a random number generator function is passed in, ensure that it only takes only argument and returns             a numpy array. The argument can be an integer or a tuple of integers. In this case, the size will be             automatically derived to save the user the trouble.          Returns         -------         simulated_data : {List[DataFrame], ndarray}             List of DataFrame with columns data containing the simulated values, volatility, containing the conditional             volatility and errors containing the errors used in the simulation          See Also         --------         UArch.simulation_horizon_required: Calculates the simulation size required Simulates data from the multiple ARCH-GARCH models.          This function is specially crafted for Monte-Carlo simulations.          Parameters         ----------         nobs: int             Length of series to simulate          reps: int             Number of repetitions          burn: int, optional             Number of values to simulate to initialize the model and remove dependence on initial values          initial_value: {ndarray, float}, optional             Either a scalar value or `max(lags)` array set of initial values to use when initializing the model.             If omitted, 0.0 is used. If array, the last column must be of the same size as the number of models          x: {ndarray, list of ndarray}, optional             If supplied as a list, this list should have the same number of elements as the number of models in the             MUArch model. Each array inside is the specified exogenous variable for that  particular model and             this must be a nobs + burn by k matrix of exogenous variables to include in the simulation. Otherwise,             leave the value as `None` to indicate no exogenous variables are used for simulation in the model.              If an array is supplied directly, it means every model has an exogenous variable associated with it. In             this case, it should be a 3 dimensional tensor where the first dimension represents the number of models.          initial_value_vol: {ndarray, float}, optional             An array or scalar to use when initializing the volatility process. If array, the last column must be of the             same size as the number of models          custom_dist: {ndarray, Callable}, optional             Optional density from which to simulate the innovations (Distribution) in the GARCH models. This is useful             when working with the copula-GARCH model where each univariate model innovations has dependence on others.             It is assumed that the values supplied are standardized [0, 1] innovations instead of the unstandardized             residuals.              The shape of the array must be at least as long as the simulation size required after accounting for burn             and type of innovation process. If unsure, use :code:`simulation_size_required` to check.              If a random number generator function is passed in, he size will be automatically derived to save the user             the trouble. However, the function must:              * take as it first argument an integer or a tuple of integer             * have other parameters that are optional             * return a numpy array          n_jobs: int or None, optional             The number of jobs to run in parallel for simulation. This is particularly useful when simulating large             number of repetitions with more than 1 dimension. None defaults to using 1 processor. Any numbers less or             equal to 0 means to use all processors. Even if a large number is used, it will be capped at the             maximum number of processors available.           Returns         -------         simulated_data : numpy array             Array containing simulated data from the Monte Carlo Simulation # setting up job processes # formulate function calls to get simulated data # switch between sequential and multiprocessing # sequential # multiprocessing Summary of fitted models          :param short: bool, default False             Whether to show short summary or full summary.         :param dp: int, default 4             Number of decimal places to show in short summary         :return: SummaryList             summary of fitted models Ensures the simulation are okay and return the formatted parameters # Format initial value # Format initial volatility value # Format exogenous regressors # convert all data frames to numpy array then move the model axis (first) to last # Format custom distribution # for monte carlo