Learn more  » Push, build, and install  RubyGems npm packages Python packages Maven artifacts PHP packages Go Modules Bower components Debian packages RPM packages NuGet packages

alkaline-ml / statsmodels   python

Repository URL to install this package:

Version: 0.11.1 

/ tsa / tests / results / results_arima.py

import os

import numpy as np
from numpy import genfromtxt
cur_dir = os.path.dirname(os.path.abspath(__file__))

path = os.path.join(cur_dir, "results_arima_forecasts.csv")
with open(path, "rb") as fd:
    forecast_results = genfromtxt(fd, names=True, delimiter=",", dtype=float)

# NOTE:
# stata gives no indication of no convergence for 112 CSS but gives a
# different answer than x12arima, gretl simply fails to converge
# redid stata with starting parameters from x12arima

# it looks like stata uses a different formula for the CSS likelihood
# they appear to be using a larger sample than R, gretl, or us.
# CSS results are therefore taken from R and gretl


class ARIMA111(object):
    def __init__(self, method="mle"):
        self.k_ar = 1
        self.k_diff = 1
        self.k_ma = 1
        if method == "mle":
            # from stata
            from .arima111_results import results

            # unpack stata results
            self.__dict__.update(results)
            self.resid = self.resid[1:]
            self.params = self.params[:-1]
            self.sigma2 = self.sigma**2
            self.aic = self.icstats[4]
            self.bic = self.icstats[5]
            self.fittedvalues = self.xb[1:]  # no idea why this initial value
            self.linear = self.y[1:]
            # stata bse are OPG
            # self.bse = np.diag(self.cov_params) ** .5

            # from gretl
            self.arroots = [1.0640 + 0j]
            self.maroots = [1.2971 + 0j]
            self.hqic = 496.8653
            self.aic_gretl = 491.5112
            self.bic_gretl = 504.7442
            self.tvalues = [4.280, 20.57, -8.590]
            self.pvalues = [1.87e-5, 5.53e-94, 8.73e-18]
            self.cov_params = [[0.0423583,   -0.00167449,    0.00262911],
                               [-0.00167449, 0.00208858,    -0.0035068],
                               [0.00262911, -0.0035068, 0.00805622]]
            self.bse = np.sqrt(np.diag(self.cov_params))
            # these bse are approx [.205811, .0457010, .0897565]

            # from stata
            # forecast = genfromtxt(open(cur_dir+"/arima111_forecasts.csv"),
            #                delimiter=",", skip_header=1, usecols=[1,2,3,4,5])
            # self.forecast = forecast[203:,1]
            # self.fcerr = forecast[203:,2]
            # self.fc_conf_int = forecast[203:,3:]

            # from gretl
            self.forecast = forecast_results['fc111c'][-25:]
            self.forecasterr = forecast_results['fc111cse'][-25:]
            self.forecast_dyn = forecast_results['fc111cdyn']
            self.forecasterr_dyn = forecast_results['fc111cdynse']
        else:
            # coefs, bse, tvalues, and pvalues taken from R because gretl
            # uses mean not constant
            self.bse = [0.21583833, 0.03844939, 0.08566390]
            self.params = [1.0087257, 0.9455393, -0.8021834]
            self.sigma2 = 0.6355913
            self.tvalues = [4.673524, 24.591788, -9.364311]
            self.pvalues = [5.464467e-06, 0, 0]
            self.cov_params = np.array([
                [0.046586183,  0.002331183, -0.004647432],
                [0.002331183,  0.001478356, -0.002726201],
                [-0.004647432, -0.002726201,  0.007338304]])

            # from gretl
            self.llf = -239.6601
            self.aic = 487.3202
            self.bic = 500.5334
            self.hqic = 492.6669
            self.arroots = [1.0578 + 0j]
            self.maroots = [1.2473 + 0j]
            # cov_params = np.array([[0.00369569, -0.00271777, 0.00269806],
            #                        [0, 0.00209573, -0.00224559],
            #                        [0, 0, 0.00342769]])
            # self.cov_params = cov_params + cov_params.T - \
            #                np.diag(np.diag(cov_params))
            # self.bse = np.sqrt(np.diag(self.cov_params))

            self.resid = [-0.015830, -0.236884, -0.093946, -0.281152,
                          -0.089983, -0.226336, -0.351666, -0.198703,
                          -0.258418, -0.259026, -0.149513, -0.325703,
                          -0.165703, -0.279229, -0.295711, -0.120018,
                          -0.289870, -0.154243, -0.348403, -0.273902,
                          -0.240894, -0.182791, -0.252930, -0.152441,
                          -0.296412, -0.128941, 0.024068, -0.243972,
                          -0.011436, -0.392437, -0.217022, -0.118190,
                          -0.133489, -0.045755, -0.169953, 0.025010,
                          -0.107754, -0.119661, 0.070794, -0.065586,
                          -0.080390, 0.007741, -0.016138, -0.235283,
                          -0.121907, -0.125546, -0.428463, -0.087713,
                          -0.298131, -0.277757, -0.261422, -0.248326,
                          -0.137826, -0.043771, 0.437100, -0.150051,
                          0.751890, 0.424180, 0.450514, 0.277089,
                          0.732583, 0.225086, -0.403648, -0.040509,
                          -0.132975, -0.112572, -0.696214, 0.003079,
                          -0.003491, -0.108758, 0.401383, -0.162302,
                          -0.141547, 0.175094, 0.245346, 0.607134, 0.519045,
                          0.248419, 0.920521, 1.097613, 0.755983, 1.271156,
                          1.216969, -0.121014, 0.340712, 0.732750, 0.068915,
                          0.603912, 0.060157, -0.803110, -1.044392, 1.040311,
                          -0.984497, -1.611668, -0.258198, -0.112970,
                          -0.091071, 0.226487, 0.097475, -0.311423, -0.061105,
                          -0.449488, 0.317277, -0.329734, -0.181248, 0.443263,
                          -2.223262, 0.096836, -0.033782, 0.456032, 0.476052,
                          0.197564, 0.263362, 0.021578, 0.216803, 0.284249,
                          0.343786, 0.196981, 0.773819, 0.169070, -0.343097,
                          0.918962, 0.096363, 0.298610, 1.571685, -0.236620,
                          -1.073822, -0.194208, -0.250742, -0.101530,
                          -0.076437, -0.056319, 0.059811, -0.041620,
                          -0.128404, -0.403446, 0.059654, -0.347208,
                          -0.095257, 0.217668, -0.015057, 0.087431, 0.275062,
                          -0.263580, -0.122746, 0.195629, 0.367272,
                          -0.184188, 0.146368, 0.127777, -0.587128,
                          -0.498538, 0.172490, -0.456741, -0.694000,
                          0.199392, -0.140634, -0.029636, 0.364818,
                          -0.097080, 0.510745, 0.230842, 0.595504, 0.709721,
                          0.012218, 0.520223, -0.445174, -0.168341,
                          -0.935465, -0.894203, 0.733417, -0.279707,
                          0.258861, 0.417969, -0.443542, -0.477955, 0.288992,
                          0.442126, 0.075826, 0.665759, 0.571509, -0.204055,
                          0.835901, -0.375693, 3.292828, -1.469299,
                          -0.122206, 0.617909, -2.250468, 0.570871, 1.166013,
                          0.079873, 0.463372, 1.981434, -0.142869, 3.023376,
                          -3.713161, -6.120150, -0.007487, 1.267027, 1.176930]

            self.linear = [
                29.3658, 29.6069, 29.6339, 29.8312, 29.8400,
                30.0663, 30.1617, 30.1187, 30.2384, 30.2990,
                30.3595, 30.5457, 30.5457, 30.7192, 30.7757,
                30.8100, 31.0399, 31.0942, 31.2984, 31.2939,
                31.3609, 31.4628, 31.6329, 31.7324, 31.9464,
                32.0089, 32.2559, 32.6940, 32.8614, 33.2924,
                33.3170, 33.5182, 33.8335, 34.1458, 34.5700,
                34.8750, 35.4078, 35.8197, 36.2292, 36.8656,
                37.3804, 37.8923, 38.5161, 39.1353, 39.5219,
                40.0255, 40.5285, 40.6877, 41.1981, 41.4778,
                41.7614, 42.0483, 42.3378, 42.7438, 43.2629,
                44.3501, 44.8481, 46.3758, 47.6495, 49.0229,
                50.2674, 52.0749, 53.4036, 54.0405, 55.0330,
                55.9126, 56.7962, 56.9969, 57.9035, 58.8088,
                59.5986, 60.9623, 61.7415, 62.5249, 63.6547,
                64.8929, 66.5810, 68.2516, 69.6795, 71.9024,
                74.4440, 76.7288, 79.6830, 82.7210, 84.3593,
                86.4672, 89.0311, 90.8961, 93.3398, 95.2031,
                96.0444, 96.4597, 99.0845, 99.5117, 99.0582,
                99.9130, 100.8911, 101.8735, 103.2025, 104.4114,
                105.1611, 106.1495, 106.6827, 108.0297, 108.6812,
                109.4567, 110.9233, 109.4032, 110.2338, 110.9440,
                112.2239, 113.6024, 114.7366, 115.9784, 116.9832,
                118.2158, 119.5562, 121.0030, 122.3262, 124.3309,
                125.7431, 126.5810, 128.8036, 130.2014, 131.8283,
                134.9366, 136.1738, 136.3942, 137.4507, 138.4015,
                139.4764, 140.5563, 141.6402, 142.8416, 143.9284,
                144.9034, 145.5403, 146.6472, 147.2953, 148.1823,
                149.4151, 150.4126, 151.5249, 152.8636, 153.6227,
                154.5044, 155.7327, 157.1842, 158.0536, 159.2722,
                160.4871, 160.8985, 161.3275, 162.4567, 162.8940,
                163.0006, 164.0406, 164.7296, 165.5352, 166.7971,
                167.5893, 169.0692, 170.3045, 171.9903, 173.8878,
                175.0798, 176.8452, 177.5683, 178.5355, 178.5942,
                178.5666, 180.2797, 180.9411, 182.1820, 183.6435,
                184.1780, 184.6110, 185.8579, 187.3242, 188.4342,
                190.2285, 192.0041, 192.9641, 195.0757, 195.9072,
                200.8693, 200.8222, 202.0821, 204.1505, 203.0031,
                204.7540, 207.2581, 208.6696, 210.5136, 214.1399,
                215.5866, 220.6022, 218.2942, 212.6785, 213.2020,
                215.2081]
            # forecasting is not any different for css
            # except you lose the first p+1 observations for in-sample
            # these results are from x-12 arima
            self.forecast = forecast_results['fc111c_css'][-25:]
            self.forecasterr = forecast_results['fc111cse_css'][-25:]
            self.forecast_dyn = forecast_results['fc111cdyn_css']
            self.forecasterr_dyn = forecast_results['fc111cdynse_css']


class ARIMA211(object):
    def __init__(self, method="mle"):
        if method == 'mle':
            # from stata
            from .arima111_results import results

            self.__dict__.update(results)
            self.resid = self.resid[1:]
            self.params = self.params[:-1]
            self.sigma2 = self.sigma**2
            self.aic = self.icstats[4]
            self.bic = self.icstats[5]
            self.fittedvalues = self.xb[1:]  # no idea why this initial value
            self.linear = self.y[1:]
            self.k_diff = 1

            # stata bse are OPG
            # self.bse = np.diag(self.cov_params) ** .5

            # from gretl
            self.arroots = [1.027 + 0j, 5.7255 + 0j]
            self.maroots = [1.1442+0j]
            self.hqic = 496.5314
            self.aic_gretl = 489.8388
            self.bic_gretl = 506.3801
            self.tvalues = [3.468, 11.14, -1.941, 12.55]
            self.pvalues = [.0005, 8.14e-29, .0522, 3.91e-36]
            cov_params = np.array([
                [0.0616906,  -0.00250187, 0.0010129,    0.00260485],
                [0, 0.0105302,   -0.00867819,   -0.00525614],
                [0, 0,         0.00759185,    0.00361962],
                [0, 0, 0,                      0.00484898]])
            self.cov_params = (
                cov_params + cov_params.T - np.diag(np.diag(cov_params)))
            self.bse = np.sqrt(np.diag(self.cov_params))
            # these bse are approx [0.248376, 0.102617, 0.0871312, 0.0696346]

            self.forecast = forecast_results['fc211c'][-25:]
            self.forecasterr = forecast_results['fc211cse'][-25:]
            self.forecast_dyn = forecast_results['fc211cdyn'][-25:]
            self.forecasterr_dyn = forecast_results['fc211cdynse'][-25:]
        else:
            from .arima211_css_results import results

            self.__dict__.update(results)
            self.resid = self.resid[1:]
            self.params = self.params[:-1]
            self.sigma2 = self.sigma**2
            self.aic = self.icstats[4]
            self.bic = self.icstats[5]
            self.fittedvalues = self.xb[1:]  # no idea why this initial value
            self.linear = self.y[1:]
            self.k_diff = 1

            # from gretl
            self.arroots = [1.0229 + 0j, 4.4501 + 0j]
            self.maroots = [1.0604 + 0j]
            self.hqic = 489.3225
            self.aic_gretl = 482.6486
            self.bic_gretl = 499.1402
            self.tvalues = [.7206, 22.54, -19.04]
            self.pvalues = [.4712, 1.52e-112, 2.19e-10, 8.00e-81]
            cov_params = np.array([
                [8.20496e-04, -0.0011992, 4.57078e-04, 0.00109907],
                [0, 0.00284432, -0.0016752, -0.00220223],
                [0, 0, 0.00119783, 0.00108868],
                [0, 0, 0, 0.00245324]])
            self.cov_params = (
                cov_params + cov_params.T - np.diag(np.diag(cov_params)))
            self.bse = np.sqrt(np.diag(self.cov_params))
            # forecasting is not any different for css
            # except you lose the first p+1 observations for in-sample
            self.forecast = forecast_results['fc111c_css'][-25:]
            self.forecasterr = forecast_results['fc111cse_css'][-25:]
            self.forecast_dyn = forecast_results['fc111cdyn_css']
            self.forecasterr_dyn = forecast_results['fc111cdynse_css']


class ARIMA112(object):
    def __init__(self, method="mle"):
        self.df_model = 3
        self.k = 5
        self.k_ar = 1
        self.k_ma = 2
        self.k_exog = 1
        self.k_diff = 1
        if method == "mle":
            from .arima112_results import results
            # from gretl
            self.arroots = [1.0324 + 0j]
            self.maroots = [1.1447 + 0j, -4.8613+0j]
            self.hqic = 495.5852
            self.aic_gretl = 488.8925
            self.bic_gretl = 505.4338
            self.tvalues = [3.454, 31.10, -7.994, -2.127]
            self.pvalues = [0.0006, 2.1e-212, 1.31e-15, .0334]
            cov_params = np.array([
                [0.0620096, -0.00172172, 0.00181301, 0.00103271],
                [0, 9.69682e-04, -9.70767e-04, -8.99814e-04],
                [0, 0, 0.00698068, -0.00443871],
                [0, 0, 0, 0.00713662]])
            self.cov_params = (
                cov_params + cov_params.T - np.diag(np.diag(cov_params)))
            self.bse = np.sqrt(np.diag(self.cov_params))

            # from gretl
            self.forecast = forecast_results['fc112c'][-25:]
            self.forecasterr = forecast_results['fc112cse'][-25:]
            self.forecast_dyn = forecast_results['fc112cdyn']
            self.forecasterr_dyn = forecast_results['fc112cdynse']

            # unpack stata results
            self.__dict__ = results
            self.resid = self.resid[1:]
            self.params = self.params[:-1]
            self.sigma2 = self.sigma**2
            self.aic = self.icstats[4]
            self.bic = self.icstats[5]
            self.fittedvalues = self.xb[1:]  # no idea why this initial value
            self.linear = self.y[1:]
            # stata bse are OPG
            # self.bse = np.diag(self.cov_params) ** .5
        else:
            # NOTE: this looks like a "hard" problem
            #  unable to replicate stata's results even with their starting
            #   values
            # unable to replicate x12 results in stata using their starting
            #   values. x-12 has better likelihood and we can replicate so
            #   use their results

            # taken from R using X12-arima values as init params
            self.bse = [0.07727588, 0.09356658, 0.10503567, 0.07727970]
            self.params = [0.9053219, -0.692412, 1.0736728, 0.1720008]
            self.sigma2 = 0.6820727
            self.tvalues = [11.715452, -7.400215, 10.221983,  2.225692]
            self.pvalues = [0, 3.791634e-12, 0, 2.716275e-02]
            self.cov_params = np.array([
                [0.0059715623, 0.001327824, -0.001592129, -0.0008061933],
                [0.0013278238, 0.008754705, -0.008024634, -0.0045933413],
                [-0.0015921293, -0.008024634,  0.011032492,  0.0072509641],
                [-0.0008061933, -0.004593341,  0.007250964,  0.0059721516]])

            # from x12arima via gretl
            # gretl did not converge for this model...
            self.llf = -246.7534
            self.nobs = 202
            # self.params = [.905322, -.692425, 1.07366, 0.172024]
            # self.sigma2 = 0.682072819129
            # self.bse = [0.0756430, 0.118440, 0.140691, 0.105266]

            self.resid = [
                -1.214477, -0.069772, -1.064510, -0.249555,
                -0.874206, -0.322177, -1.003579, -0.310040, -0.890506,
Loading ...