https://stackoverflow.com/questions/44301263/python-arma-mle-implementing-algorithms-from-literature
已知外國有人他自己看 literature 實做出來 ARIMA ，這樣就不用 import 別人早就寫好的 python statsmodel 的 ARIMA 類別與方法，讓自己讀懂 literature ，知道 ARIMA 究竟想要做什麼，實現自己想要的情境，請問下述的 code 具體執行了什麼？謝謝！

import numpy as np
import statsmodels.api as sm
import statsmodels.tsa.arima_model
import scipy.optimize

class ARMA(object):
    def __init__(self, endo, nar, nma):
        np.random.seed(0)
        # endogenous variables
        self.endo = endo
        # Number of AR terms
        self.nar = nar
        # Number of MA terms
        self.nma = nma
        # "Dimension" of the ARMA fit
        self.dim = max(nar, nma+1)
        # Current ARMA parameters
        self.params = np.zeros(self.nar+self.nma, dtype='float')

    def __g(self, ma_params):
        '''
        Build MA parameter vector
        '''
        g = np.zeros(self.dim, dtype='float')
        g[0] = 1.0
        if self.nma > 0:
            g[1:self.nma+1] = ma_params
        return g

    def __F(self, ar_params):
        '''
        Build AR parameter matrix
        '''
        F = np.zeros((self.dim, self.dim), dtype='float')
        F[:self.nar, 0] =  ar_params
        for i in xrange(1, self.dim):
            F[i-1, i] = 1.0
        return F

    def __initial_P(self, R, T):
        '''
        Solve for initial P matrix

        Solves P = TPT' + RR'
        '''
        v = np.zeros(self.dim*self.dim, dtype='float')
        for i in xrange(self.dim):
            for j in xrange(self.dim):
                v[i+j*self.dim] = R[i]*R[j]
        R = np.array([R])
        S = np.identity(self.dim**2, dtype='float')-np.kron(T, T)
        V = np.outer(R, R).ravel('F')
        Pmat = np.linalg.solve(S,V).reshape(self.dim, self.dim, order='F')
        return Pmat

    def __likelihood(self, params):
        '''
        Compute log likehood for a parameter vector

        Implements the Pearlman 1980 algorithm
        '''
        # these checks are pilfered from statsmodels
        if self.nar > 0 and not np.all(np.abs(np.roots(np.r_[1, -params[:self.nar]]))<1):
            print 'AR coefficients are not stationary'
        if self.nma > 0 and not np.all(np.abs(np.roots(np.r_[1, -params[-self.nma:]]))<1):
           print 'MA coefficients are not stationary'
        ar_params = params[:self.nar]
        ma_params = params[-self.nma:]
        g = self.__g(ma_params)
        F = self.__F(ar_params)
        w = self.endo
        P = self.__initial_P(g, F)
        n = len(w)
        z = np.zeros(self.dim, dtype='float')
        R = np.zeros(n, dtype='float')
        a = np.zeros(n, dtype='float')
        K = np.dot(F, P[:, 0])
        L = K.copy()
        R[0] = P[0, 0]
        for i in xrange(1, n):
            a[i-1] = w[i-1] - z[0]
            z = np.dot(F, z) + K*(a[i-1]/R[i-1])
            Kupdate = -(L[0]/R[i-1])*np.dot(F, L)
            Rupdate = -L[0]*L[0]/R[i-1]
            P -= np.outer(L, L)/R[i-1]
            L = np.dot(F, L) - (L[0]/R[i-1])*K
            K += Kupdate
            R[i] = R[i-1] + Rupdate
            if np.abs(R[i] - 1.0) < 1e-9:
                R[i:] = 1.0
                break
        for j in xrange(i, n):
            a[j] = w[j] - z[0]
            z = np.dot(F, z) + K*(a[i-1]/R[i-1])
        likelihood = 0.0
        for i in xrange(n):
            likelihood += np.log(R[i])
        likelihood *= -0.5
        ssum = 0.0
        for i in xrange(n):
            ssum += a[i]*a[i]/R[i]
        likelihood += -0.5*n*np.log(ssum)
        return likelihood

    def fit(self):
        '''
        Fit the ARMA model by minimizing the loglikehood

        Uses scipy.optimize.minimize
        '''
        sm_arma = statsmodels.tsa.arima_model.ARMA(endog=self.endo, order=(self.nar, self.nma, 0))
        params = statsmodels.tsa.arima_model.ARMA._fit_start_params_hr(sm_arma, order=(self.nar, self.nma, 0))
        opt = scipy.optimize.minimize(fun=self.__likelihood, x0=params, method='L-BFGS-B')
        print opt

# Test the code on statsmodels sunspots data
nar = 2
nma = 1

endo = sm.datasets.sunspots.load_pandas().data['SUNACTIVITY'].tolist()
arma = ARMA(endo=endo, nar=nar, nma=nma)
arma.fit()