<s> '''<STR_LIT>''' <EOL> from __future__ import print_function <EOL> from statsmodels . compat . python import range <EOL> import numpy as np <EOL> from scipy import signal , optimize , linalg <EOL> def arma_generate_sample ( ar , ma , nsample , sigma = <NUM_LIT:1> , distrvs = np . random . randn , <EOL> burnin = <NUM_LIT:0> ) : <EOL> """<STR_LIT>""" <EOL> eta = sigma * distrvs ( nsample + burnin ) <EOL> return signal . lfilter ( ma , ar , eta ) [ burnin : ] <EOL> def arma_acovf ( ar , ma , nobs = <NUM_LIT:10> ) : <EOL> '''<STR_LIT>''' <EOL> if np . abs ( np . sum ( ar ) - <NUM_LIT:1> ) > <NUM_LIT> : <EOL> nobs_ir = max ( <NUM_LIT:1000> , <NUM_LIT:2> * nobs ) <EOL> else : <EOL> nobs_ir = max ( <NUM_LIT:100> , <NUM_LIT:2> * nobs ) <EOL> ir = arma_impulse_response ( ar , ma , nobs = nobs_ir ) <EOL> while ir [ - <NUM_LIT:1> ] > <NUM_LIT:5> * <NUM_LIT> : <EOL> nobs_ir *= <NUM_LIT:10> <EOL> ir = arma_impulse_response ( ar , ma , nobs = nobs_ir ) <EOL> if nobs_ir > <NUM_LIT> and nobs < <NUM_LIT> : <EOL> acovf = np . array ( [ np . dot ( ir [ : nobs - t ] , ir [ t : nobs ] ) <EOL> for t in range ( nobs ) ] ) <EOL> else : <EOL> acovf = np . correlate ( ir , ir , '<STR_LIT>' ) [ len ( ir ) - <NUM_LIT:1> : ] <EOL> return acovf [ : nobs ] <EOL> def arma_acf ( ar , ma , nobs = <NUM_LIT:10> ) : <EOL> '''<STR_LIT>''' <EOL> acovf = arma_acovf ( ar , ma , nobs ) <EOL> return acovf / acovf [ <NUM_LIT:0> ] <EOL> def arma_pacf ( ar , ma , nobs = <NUM_LIT:10> ) : <EOL> '''<STR_LIT>''' <EOL> apacf = np . zeros ( nobs ) <EOL> acov = arma_acf ( ar , ma , nobs = nobs + <NUM_LIT:1> ) <EOL> apacf [ <NUM_LIT:0> ] = <NUM_LIT:1.> <EOL> for k in range ( <NUM_LIT:2> , nobs + <NUM_LIT:1> ) : <EOL> r = acov [ : k ] <EOL> apacf [ k - <NUM_LIT:1> ] = linalg . solve ( linalg . toeplitz ( r [ : - <NUM_LIT:1> ] ) , r [ <NUM_LIT:1> : ] ) [ - <NUM_LIT:1> ] <EOL> return apacf <EOL> def arma_periodogram ( ar , ma , worN = None , whole = <NUM_LIT:0> ) : <EOL> '''<STR_LIT>''' <EOL> w , h = signal . freqz ( ma , ar , worN = worN , whole = whole ) <EOL> sd = np . abs ( h ) ** <NUM_LIT:2> / np . sqrt ( <NUM_LIT:2> * np . pi ) <EOL> if np . sum ( np . isnan ( h ) ) > <NUM_LIT:0> : <EOL> print ( '<STR_LIT>' ) <EOL> return w , sd <EOL> def arma_impulse_response ( ar , ma , nobs = <NUM_LIT:100> ) : <EOL> '''<STR_LIT>''' <EOL> impulse = np . zeros ( nobs ) <EOL> impulse [ <NUM_LIT:0> ] = <NUM_LIT:1.> <EOL> return signal . lfilter ( ma , ar , impulse ) <EOL> arma2ma = arma_impulse_response <EOL> def arma2ar ( ar , ma , nobs = <NUM_LIT:100> ) : <EOL> '''<STR_LIT>''' <EOL> return arma_impulse_response ( ma , ar , nobs = nobs ) <EOL> def ar2arma ( ar_des , p , q , n = <NUM_LIT:20> , mse = '<STR_LIT>' , start = None ) : <EOL> '''<STR_LIT>''' <EOL> def msear_err ( arma , ar_des ) : <EOL> ar , ma = np . r_ [ <NUM_LIT:1> , arma [ : p - <NUM_LIT:1> ] ] , np . r_ [ <NUM_LIT:1> , arma [ p - <NUM_LIT:1> : ] ] <EOL> ar_approx = arma_impulse_response ( ma , ar , n ) <EOL> return ( ar_des - ar_approx ) <EOL> if start is None : <EOL> arma0 = np . r_ [ - <NUM_LIT> * np . ones ( p - <NUM_LIT:1> ) , np . zeros ( q - <NUM_LIT:1> ) ] <EOL> else : <EOL> arma0 = start <EOL> res = optimize . leastsq ( msear_err , arma0 , ar_des , maxfev = <NUM_LIT> ) <EOL> arma_app = np . atleast_1d ( res [ <NUM_LIT:0> ] ) <EOL> ar_app = np . r_ [ <NUM_LIT:1> , arma_app [ : p - <NUM_LIT:1> ] ] , <EOL> ma_app = np . r_ [ <NUM_LIT:1> , arma_app [ p - <NUM_LIT:1> : ] ] <EOL> return ar_app , ma_app , res <EOL> def lpol2index ( ar ) : <EOL> '''<STR_LIT>''' <EOL> ar = np . asarray ( ar ) <EOL> index = np . nonzero ( ar ) [ <NUM_LIT:0> ] <EOL> coeffs = ar [ index ] <EOL> return coeffs , index <EOL> def index2lpol ( coeffs , index ) : <EOL> '''<STR_LIT>''' <EOL> n = max ( index ) <EOL> ar = np . zeros ( n ) <EOL> ar [ index ] = coeffs <EOL> return ar <EOL> def lpol_fima ( d , n = <NUM_LIT:20> ) : <EOL> '''<STR_LIT>''' <EOL> from scipy . special import gammaln <EOL> j = np . arange ( n ) <EOL> return np . exp ( gammaln ( d + j ) - gammaln ( j + <NUM_LIT:1> ) - gammaln ( d ) ) <EOL> def lpol_fiar ( d , n = <NUM_LIT:20> ) : <EOL> '''<STR_LIT>''' <EOL> from scipy . special import gammaln <EOL> j = np . arange ( n ) <EOL> ar = - np . exp ( gammaln ( - d + j ) - gammaln ( j + <NUM_LIT:1> ) - gammaln ( - d ) ) <EOL> ar [ <NUM_LIT:0> ] = <NUM_LIT:1> <EOL> return ar <EOL> def lpol_sdiff ( s ) : <EOL> '''<STR_LIT>''' <EOL> return [ <NUM_LIT:1> ] + [ <NUM_LIT:0> ] * ( s - <NUM_LIT:1> ) + [ - <NUM_LIT:1> ] <EOL> def deconvolve ( num , den , n = None ) : <EOL> """<STR_LIT>""" <EOL> num = np . atleast_1d ( num ) <EOL> den = np . atleast_1d ( den ) <EOL> N = len ( num ) <EOL> D = len ( den ) <EOL> if D > N and n is None : <EOL> quot = [ ] <EOL> rem = num <EOL> else : <EOL> if n is None : <EOL> n = N - D + <NUM_LIT:1> <EOL> input = np . zeros ( n , float ) <EOL> input [ <NUM_LIT:0> ] = <NUM_LIT:1> <EOL> quot = signal . lfilter ( num , den , input ) <EOL> num_approx = signal . convolve ( den , quot , mode = '<STR_LIT>' ) <EOL> if len ( num ) < len ( num_approx ) : <EOL> num = np . concatenate ( ( num , np . zeros ( len ( num_approx ) - len ( num ) ) ) ) <EOL> rem = num - num_approx <EOL> return quot , rem <EOL> class ArmaProcess ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , ar , ma , nobs = <NUM_LIT:100> ) : <EOL> self . ar = np . asarray ( ar ) <EOL> self . ma = np . asarray ( ma ) <EOL> self . arcoefs = - self . ar [ <NUM_LIT:1> : ] <EOL> self . macoefs = self . ma [ <NUM_LIT:1> : ] <EOL> self . arpoly = np . polynomial . Polynomial ( self . ar ) <EOL> self . mapoly = np . polynomial . Polynomial ( self . ma ) <EOL> self . nobs = nobs <EOL> @ classmethod <EOL> def from_coeffs ( cls , arcoefs , macoefs , nobs = <NUM_LIT:100> ) : <EOL> """<STR_LIT>""" <EOL> return cls ( np . r_ [ <NUM_LIT:1> , - arcoefs ] , np . r_ [ <NUM_LIT:1> , macoefs ] , nobs = nobs ) <EOL> @ classmethod <EOL> def from_estimation ( cls , model_results , nobs = None ) : <EOL> """<STR_LIT>""" <EOL> arcoefs = model_results . arparams <EOL> macoefs = model_results . maparams <EOL> nobs = nobs or model_results . nobs <EOL> return cls ( np . r_ [ <NUM_LIT:1> , - arcoefs ] , np . r_ [ <NUM_LIT:1> , macoefs ] , nobs = nobs ) <EOL> def __mul__ ( self , oth ) : <EOL> if isinstance ( oth , self . __class__ ) : <EOL> ar = ( self . arpoly * oth . arpoly ) . coef <EOL> ma = ( self . mapoly * oth . mapoly ) . coef <EOL> else : <EOL> try : <EOL> aroth , maoth = oth <EOL> arpolyoth = np . polynomial . Polynomial ( aroth ) <EOL> mapolyoth = np . polynomial . Polynomial ( maoth ) <EOL> ar = ( self . arpoly * arpolyoth ) . coef <EOL> ma = ( self . mapoly * mapolyoth ) . coef <EOL> except : <EOL> print ( '<STR_LIT>' ) <EOL> raise <EOL> return self . __class__ ( ar , ma , nobs = self . nobs ) <EOL> def __repr__ ( self ) : <EOL> return '<STR_LIT>' % ( self . ar . tolist ( ) , <EOL> self . ma . tolist ( ) , <EOL> self . nobs ) <EOL> def __str__ ( self ) : <EOL> return '<STR_LIT>' % ( self . ar . tolist ( ) , <EOL> self . ma . tolist ( ) ) <EOL> def acovf ( self , nobs = None ) : <EOL> nobs = nobs or self . nobs <EOL> return arma_acovf ( self . ar , self . ma , nobs = nobs ) <EOL> acovf . __doc__ = arma_acovf . __doc__ <EOL> def acf ( self , nobs = None ) : <EOL> nobs = nobs or self . nobs <EOL> return arma_acf ( self . ar , self . ma , nobs = nobs ) <EOL> acf . __doc__ = arma_acf . __doc__ <EOL> def pacf ( self , nobs = None ) : <EOL> nobs = nobs or self . nobs <EOL> return arma_pacf ( self . ar , self . ma , nobs = nobs ) <EOL> pacf . __doc__ = arma_pacf . __doc__ <EOL> def periodogram ( self , nobs = None ) : <EOL> nobs = nobs or self . nobs <EOL> return arma_periodogram ( self . ar , self . ma , worN = nobs ) <EOL> periodogram . __doc__ = arma_periodogram . __doc__ <EOL> def impulse_response ( self , nobs = None ) : <EOL> nobs = nobs or self . nobs <EOL> return arma_impulse_response ( self . ar , self . ma , worN = nobs ) <EOL> impulse_response . __doc__ = arma_impulse_response . __doc__ <EOL> def arma2ma ( self , nobs = None ) : <EOL> nobs = nobs or self . nobs <EOL> return arma2ma ( self . ar , self . ma , nobs = nobs ) <EOL> arma2ma . __doc__ = arma2ma . __doc__ <EOL> def arma2ar ( self , nobs = None ) : <EOL> nobs = nobs or self . nobs <EOL> return arma2ar ( self . ar , self . ma , nobs = nobs ) <EOL> arma2ar . __doc__ = arma2ar . __doc__ <EOL> @ property <EOL> def arroots ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . arpoly . roots ( ) <EOL> @ property <EOL> def maroots ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . mapoly . roots ( ) <EOL> @ property <EOL> def isstationary ( self ) : <EOL> '''<STR_LIT>''' <EOL> if np . all ( np . abs ( self . arroots ) > <NUM_LIT:1> ) : <EOL> return True <EOL> else : <EOL> return False <EOL> @ property <EOL> def isinvertible ( self ) : <EOL> '''<STR_LIT>''' <EOL> if np . all ( np . abs ( self . maroots ) > <NUM_LIT:1> ) : <EOL> return True <EOL> else : <EOL> return False <EOL> def invertroots ( self , retnew = False ) : <EOL> '''<STR_LIT>''' <EOL> pr = self . ma_roots ( ) <EOL> insideroots = np . abs ( pr ) < <NUM_LIT:1> <EOL> if insideroots . any ( ) : <EOL> pr [ np . abs ( pr ) < <NUM_LIT:1> ] = <NUM_LIT:1.> / pr [ np . abs ( pr ) < <NUM_LIT:1> ] <EOL> pnew = np . polynomial . Polynomial . fromroots ( pr ) <EOL> mainv = pnew . coef / pnew . coef [ <NUM_LIT:0> ] <EOL> wasinvertible = False <EOL> else : <EOL> mainv = self . ma <EOL> wasinvertible = True <EOL> if retnew : <EOL> return self . __class__ ( self . ar , mainv , nobs = self . nobs ) <EOL> else : <EOL> return mainv , wasinvertible <EOL> def generate_sample ( self , nsample = <NUM_LIT:100> , scale = <NUM_LIT:1.> , distrvs = None , axis = <NUM_LIT:0> , <EOL> burnin = <NUM_LIT:0> ) : <EOL> '''<STR_LIT>''' <EOL> if distrvs is None : <EOL> distrvs = np . random . normal <EOL> if np . ndim ( nsample ) == <NUM_LIT:0> : <EOL> nsample = [ nsample ] <EOL> if burnin : <EOL> newsize = list ( nsample ) <EOL> newsize [ axis ] += burnin <EOL> newsize = tuple ( newsize ) <EOL> fslice = [ slice ( None ) ] * len ( newsize ) <EOL> fslice [ axis ] = slice ( burnin , None , None ) <EOL> fslice = tuple ( fslice ) <EOL> else : <EOL> newsize = tuple ( nsample ) <EOL> fslice = tuple ( [ slice ( None ) ] * np . ndim ( newsize ) ) <EOL> eta = scale * distrvs ( size = newsize ) <EOL> return signal . lfilter ( self . ma , self . ar , eta , axis = axis ) [ fslice ] <EOL> __all__ = [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , <EOL> '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , <EOL> '<STR_LIT>' , '<STR_LIT>' ] <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> ar = [ <NUM_LIT:1> , - <NUM_LIT> ] <EOL> ma = [ <NUM_LIT:1.0> ] <EOL> eta = <NUM_LIT:0.1> * np . random . randn ( <NUM_LIT:1000> ) <EOL> yar1 = signal . lfilter ( ar , ma , eta ) <EOL> print ( "<STR_LIT>" ) <EOL> arest = ARIMAProcess ( yar1 ) <EOL> rhohat , cov_x , infodict , mesg , ier = arest . fit ( ( <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:1> ) ) <EOL> print ( rhohat ) <EOL> print ( cov_x ) <EOL> print ( "<STR_LIT>" ) <EOL> ar = [ <NUM_LIT:1.0> , - <NUM_LIT> ] <EOL> ma = [ <NUM_LIT:1.0> , <NUM_LIT:0.5> ] <EOL> y1 = arest . generate_sample ( ar , ma , <NUM_LIT:1000> , <NUM_LIT:0.1> ) <EOL> arest = ARIMAProcess ( y1 ) <EOL> rhohat1 , cov_x1 , infodict , mesg , ier = arest . fit ( ( <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:1> ) ) <EOL> print ( rhohat1 ) <EOL> print ( cov_x1 ) <EOL> err1 = arest . errfn ( x = y1 ) <EOL> print ( np . var ( err1 ) ) <EOL> import statsmodels . api as sm <EOL> print ( sm . regression . yule_walker ( y1 , order = <NUM_LIT:2> , inv = True ) ) <EOL> print ( "<STR_LIT>" ) <EOL> nsample = <NUM_LIT:1000> <EOL> ar = [ <NUM_LIT:1.0> , - <NUM_LIT> , - <NUM_LIT:0.1> ] <EOL> ma = [ <NUM_LIT:1.0> , <NUM_LIT> , <NUM_LIT> ] <EOL> y2 = ARIMA . generate_sample ( ar , ma , nsample , <NUM_LIT:0.1> ) <EOL> arest2 = ARIMAProcess ( y2 ) <EOL> rhohat2 , cov_x2 , infodict , mesg , ier = arest2 . fit ( ( <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:2> ) ) <EOL> print ( rhohat2 ) <EOL> print ( cov_x2 ) <EOL> err2 = arest . errfn ( x = y2 ) <EOL> print ( np . var ( err2 ) ) <EOL> print ( arest2 . rhoy ) <EOL> print ( arest2 . rhoe ) <EOL> print ( "<STR_LIT:true>" ) <EOL> print ( ar ) <EOL> print ( ma ) <EOL> rhohat2a , cov_x2a , infodict , mesg , ier = arest2 . fit ( ( <NUM_LIT:2> , <NUM_LIT:0> , <NUM_LIT:2> ) ) <EOL> print ( rhohat2a ) <EOL> print ( cov_x2a ) <EOL> err2a = arest . errfn ( x = y2 ) <EOL> print ( np . var ( err2a ) ) <EOL> print ( arest2 . rhoy ) <EOL> print ( arest2 . rhoe ) <EOL> print ( "<STR_LIT:true>" ) <EOL> print ( ar ) <EOL> print ( ma ) <EOL> print ( sm . regression . yule_walker ( y2 , order = <NUM_LIT:2> , inv = True ) ) <EOL> print ( "<STR_LIT>" ) <EOL> nsample = <NUM_LIT:1000> <EOL> ar = [ <NUM_LIT:1.0> ] <EOL> ma = [ <NUM_LIT:1.0> , <NUM_LIT:0.5> , <NUM_LIT> ] <EOL> y3 = ARIMA . generate_sample ( ar , ma , nsample , <NUM_LIT> ) <EOL> arest20 = ARIMAProcess ( y3 ) <EOL> rhohat3 , cov_x3 , infodict , mesg , ier = arest20 . fit ( ( <NUM_LIT:2> , <NUM_LIT:0> , <NUM_LIT:0> ) ) <EOL> print ( rhohat3 ) <EOL> print ( cov_x3 ) <EOL> err3 = arest20 . errfn ( x = y3 ) <EOL> print ( np . var ( err3 ) ) <EOL> print ( np . sqrt ( np . dot ( err3 , err3 ) / nsample ) ) <EOL> print ( arest20 . rhoy ) <EOL> print ( arest20 . rhoe ) <EOL> print ( "<STR_LIT:true>" ) <EOL> print ( ar ) <EOL> print ( ma ) <EOL> rhohat3a , cov_x3a , infodict , mesg , ier = arest20 . fit ( ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:2> ) ) <EOL> print ( rhohat3a ) <EOL> print ( cov_x3a ) <EOL> err3a = arest20 . errfn ( x = y3 ) <EOL> print ( np . var ( err3a ) ) <EOL> print ( np . sqrt ( np . dot ( err3a , err3a ) / nsample ) ) <EOL> print ( arest20 . rhoy ) <EOL> print ( arest20 . rhoe ) <EOL> print ( "<STR_LIT:true>" ) <EOL> print ( ar ) <EOL> print ( ma ) <EOL> print ( sm . regression . yule_walker ( y3 , order = <NUM_LIT:2> , inv = True ) ) <EOL> print ( "<STR_LIT>" ) <EOL> nsample = <NUM_LIT:1000> <EOL> ar = [ <NUM_LIT:1.0> , - <NUM_LIT> , <NUM_LIT> ] <EOL> ma = [ <NUM_LIT:1.0> ] <EOL> y4 = ARIMA . generate_sample ( ar , ma , nsample ) <EOL> arest02 = ARIMAProcess ( y4 ) <EOL> rhohat4 , cov_x4 , infodict , mesg , ier = arest02 . fit ( ( <NUM_LIT:2> , <NUM_LIT:0> , <NUM_LIT:0> ) ) <EOL> print ( rhohat4 ) <EOL> print ( cov_x4 ) <EOL> err4 = arest02 . errfn ( x = y4 ) <EOL> print ( np . var ( err4 ) ) <EOL> sige = np . sqrt ( np . dot ( err4 , err4 ) / nsample ) <EOL> print ( sige ) <EOL> print ( sige * np . sqrt ( np . diag ( cov_x4 ) ) ) <EOL> print ( np . sqrt ( np . diag ( cov_x4 ) ) ) <EOL> print ( arest02 . rhoy ) <EOL> print ( arest02 . rhoe ) <EOL> print ( "<STR_LIT:true>" ) <EOL> print ( ar ) <EOL> print ( ma ) <EOL> rhohat4a , cov_x4a , infodict , mesg , ier = arest02 . fit ( ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:2> ) ) <EOL> print ( rhohat4a ) <EOL> print ( cov_x4a ) <EOL> err4a = arest02 . errfn ( x = y4 ) <EOL> print ( np . var ( err4a ) ) <EOL> sige = np . sqrt ( np . dot ( err4a , err4a ) / nsample ) <EOL> print ( sige ) <EOL> print ( sige * np . sqrt ( np . diag ( cov_x4a ) ) ) <EOL> print ( np . sqrt ( np . diag ( cov_x4a ) ) ) <EOL> print ( arest02 . rhoy ) <EOL> print ( arest02 . rhoe ) <EOL> print ( "<STR_LIT:true>" ) <EOL> print ( ar ) <EOL> print ( ma ) <EOL> import statsmodels . api as sm <EOL> print ( sm . regression . yule_walker ( y4 , order = <NUM_LIT:2> , method = '<STR_LIT>' , inv = True ) ) <EOL> import matplotlib . pyplot as plt <EOL> plt . plot ( arest2 . forecast ( ) [ - <NUM_LIT:100> : ] ) <EOL> ar1 , ar2 = ( [ <NUM_LIT:1> , - <NUM_LIT> ] , [ <NUM_LIT:1> , <NUM_LIT:0.5> ] ) <EOL> ar2 = [ <NUM_LIT:1> , - <NUM_LIT:1> ] <EOL> lagpolyproduct = np . convolve ( ar1 , ar2 ) <EOL> print ( deconvolve ( lagpolyproduct , ar2 , n = None ) ) <EOL> print ( signal . deconvolve ( lagpolyproduct , ar2 ) ) <EOL> print ( deconvolve ( lagpolyproduct , ar2 , n = <NUM_LIT:10> ) ) </s>
