arfit.py
def arord(R, m, mcor, ne, pmin, pmax):
    '''Evaluates criteria for selecting the order of an AR model.
    It returns approximate values of the order selection criteria 
        SBC and FPE for models of order pmin:pmax.
    @R: An upper triangular matrix
    @m: dimension of state vectors
    @mcor: a flag indicates whether or not an intercept vector
        is being fitted
    @ne: number of block equations of size m
    @pmin: minimum order
    @pmax: maximum order
    @return: SBC and FPE estimation results, in vectors
    '''
    imax = pmax - pmin + 1 # maximum index of output vectors

    sbc = np.zeros((imax,)).T
    fpe = np.zeros((imax,)).T
    logdp = np.zeros((imax,)).T
    np_vec = np.zeros((imax,), dtype=np.int).T
    np_vec[imax-1] = int(m * pmax + mcor)

    # Get lower right triangle R22 of R: 
    #     | R11  R12 |
    # R=  |          |
    #     | 0    R22 |
    R22 = R[np_vec[imax-1]-1:np_vec[imax-1]+m-1, np_vec[imax-1]-1:np_vec[imax-1]+m-1]
    # get inverse of residual cross-product matrix for model
    # of order pmax
    invR22 = numpy.linalg.inv(R22)
    Mp = invR22 * invR22.T

    # For order selection, get determinant of residual cross-product matrix
    #   logdp = log det(residual cross-product matrix)
    logdp[imax-1] = 2 * np.log(np.abs(R22))

    # Compute approximate order selection criteria for models of 
    #   order pmin:pmax
    i = imax - 1
    for p in reversed(range(pmin, pmax)):
        np_vec[i-1] = m * p + mcor
        if p < pmax:
            # Downdate determinant of residual cross-product matrix
            # Rp: Part of R to be added to Cholesky factor of covariance matrix
            Rp = R[np_vec[i-1]:np_vec[i-1]+m, np_vec[imax-1]:np_vec[imax-1]+m]

            # Get Mp, the downdated inverse of the residual cross-product
            # matrix, using the Woodbury formula
            L = numpy.linalg.cholesky(np.eye(m) + Rp * Mp * Rp.T)
            # note: here since we only deal with 1-dimensional data,
            # in order to solve LN = Rp*Mp (where in matlab we call 'L \ Rp*Mp'),
            # we simply perform division to get value of N
            N = Rp * Mp / L
            Mp -= N**2

            # get downgraded logarithm of determinant
            logdp[i-1] = logdp[i] + 2 * np.log(np.abs(L))
    
        # SBC
        sbc[i-1] = logdp[i-1] / m - np.log(ne) * (ne-np_vec[i-1]) / ne

        # FPE
        fpe[i-1] = logdp[i-1] / m - np.log(ne * (ne-np_vec[i-1]) / (ne+np_vec[i-1]))
        
        i -= 1
    
    return sbc, fpe
Evaluates criteria for selecting the order of an AR model.
    It returns approximate values of the order selection criteria 
        SBC and FPE for models of order pmin:pmax.
    @R: An upper triangular matrix
    @m: dimension of state vectors
    @mcor: a flag indicates whether or not an intercept vector
        is being fitted
    @ne: number of block equations of size m
    @pmin: minimum order
    @pmax: maximum order
    @return: SBC and FPE estimation results, in vectors
    
Evaluates criteria for selecting the order of an AR model.
It returns approximate values of the order selection criteria
SBC and FPE for models of order pmin:pmax.
Evaluates criteria for selecting the order of an AR model . It returns approximate values of the order selection criteria SBC and FPE for models of order pmin : pmax .
def arord(R, m, mcor, ne, pmin, pmax):
    imax = pmax - pmin + 1 
    sbc = np.zeros((imax,)).T
    fpe = np.zeros((imax,)).T
    logdp = np.zeros((imax,)).T
    np_vec = np.zeros((imax,), dtype=np.int).T
    np_vec[imax-1] = int(m * pmax + mcor)
    R22 = R[np_vec[imax-1]-1:np_vec[imax-1]+m-1, np_vec[imax-1]-1:np_vec[imax-1]+m-1]
    invR22 = numpy.linalg.inv(R22)
    Mp = invR22 * invR22.T
    logdp[imax-1] = 2 * np.log(np.abs(R22))
    i = imax - 1
    for p in reversed(range(pmin, pmax)):
        np_vec[i-1] = m * p + mcor
        if p < pmax:
            Rp = R[np_vec[i-1]:np_vec[i-1]+m, np_vec[imax-1]:np_vec[imax-1]+m]
            L = numpy.linalg.cholesky(np.eye(m) + Rp * Mp * Rp.T)
            N = Rp * Mp / L
            Mp -= N**2
            logdp[i-1] = logdp[i] + 2 * np.log(np.abs(L))
        sbc[i-1] = logdp[i-1] / m - np.log(ne) * (ne-np_vec[i-1]) / ne
        fpe[i-1] = logdp[i-1] / m - np.log(ne * (ne-np_vec[i-1]) / (ne+np_vec[i-1]))
        i -= 1
    return sbc, fpe
def arord ( R , m , mcor , ne , pmin , pmax ) : imax = pmax - pmin + 1 sbc = np . zeros ( ( imax , ) ) . T fpe = np . zeros ( ( imax , ) ) . T logdp = np . zeros ( ( imax , ) ) . T np_vec = np . zeros ( ( imax , ) , dtype = np . int ) . T np_vec [ imax - 1 ] = int ( m * pmax + mcor ) R22 = R [ np_vec [ imax - 1 ] - 1 : np_vec [ imax - 1 ] + m - 1 , np_vec [ imax - 1 ] - 1 : np_vec [ imax - 1 ] + m - 1 ] invR22 = numpy . linalg . inv ( R22 ) Mp = invR22 * invR22 . T logdp [ imax - 1 ] = 2 * np . log ( np . abs ( R22 ) ) i = imax - 1 for p in reversed ( range ( pmin , pmax ) ) : np_vec [ i - 1 ] = m * p + mcor if p < pmax : Rp = R [ np_vec [ i - 1 ] : np_vec [ i - 1 ] + m , np_vec [ imax - 1 ] : np_vec [ imax - 1 ] + m ] L = numpy . linalg . cholesky ( np . eye ( m ) + Rp * Mp * Rp . T ) N = Rp * Mp / L Mp -= N ** 2 logdp [ i - 1 ] = logdp [ i ] + 2 * np . log ( np . abs ( L ) ) sbc [ i - 1 ] = logdp [ i - 1 ] / m - np . log ( ne ) * ( ne - np_vec [ i - 1 ] ) / ne fpe [ i - 1 ] = logdp [ i - 1 ] / m - np . log ( ne * ( ne - np_vec [ i - 1 ] ) / ( ne + np_vec [ i - 1 ] ) ) i -= 1 return sbc , fpe
Evaluates criteria for selecting the order of an AR model.
Evaluates criteria for selecting the order of an AR model .
'''Evaluates criteria for selecting the order of an AR model.
    It returns approximate values of the order selection criteria 
        SBC and FPE for models of order pmin:pmax.
    @R: An upper triangular matrix
    @m: dimension of state vectors
    @mcor: a flag indicates whether or not an intercept vector
        is being fitted
    @ne: number of block equations of size m
    @pmin: minimum order
    @pmax: maximum order
    @return: SBC and FPE estimation results, in vectors
    ''' # maximum index of output vectors # Get lower right triangle R22 of R:  #     | R11  R12 | # R=  |          | #     | 0    R22 | # get inverse of residual cross-product matrix for model # of order pmax # For order selection, get determinant of residual cross-product matrix #   logdp = log det(residual cross-product matrix) # Compute approximate order selection criteria for models of  #   order pmin:pmax # Downdate determinant of residual cross-product matrix # Rp: Part of R to be added to Cholesky factor of covariance matrix # Get Mp, the downdated inverse of the residual cross-product # matrix, using the Woodbury formula # note: here since we only deal with 1-dimensional data, # in order to solve LN = Rp*Mp (where in matlab we call 'L \ Rp*Mp'), # we simply perform division to get value of N # get downgraded logarithm of determinant # SBC # FPE
R
m
mcor
ne
pmin
pmax