0
# ARMA Methods
1

2
Auto-regressive moving average (ARMA) parameter estimation and power spectral density computation for combined AR-MA models and pure MA models. These methods extend AR modeling by including a moving average component, providing more flexible spectral modeling capabilities.
3

4
## Capabilities
5

6
### ARMA Parameter Estimation
7

8
Estimates both AR and MA parameters for ARMA(P,Q) models using modified Yule-Walker technique.
9

10
```python { .api }
11
def arma_estimate(X, P, Q, lag):
12
    """
13
    Estimate ARMA(P,Q) parameters using modified Yule-Walker technique.
14
    
15
    Parameters:
16
    - X: array-like, input time series data
17
    - P: int, AR model order
18
    - Q: int, MA model order  
19
    - lag: int, maximum lag for correlation computation
20
    
21
    Returns:
22
    tuple: (AR parameters array, MA parameters array, white noise variance float)
23
    """
24
```
25

26
### Moving Average Estimation
27

28
Estimates MA parameters using long AR model approximation approach.
29

30
```python { .api }
31
def ma(X, Q, M):
32
    """
33
    Moving average estimator using long AR model approach.
34
    
35
    Parameters:
36
    - X: array-like, input time series data
37
    - Q: int, MA model order
38
    - M: int, intermediate AR model order (should be >> Q)
39
    
40
    Returns:
41
    tuple: (MA parameters array, white noise variance estimate float)
42
    """
43
```
44

45
### ARMA to PSD Conversion
46

47
Computes power spectral density given ARMA model parameters, supporting pure AR, pure MA, or combined ARMA models.
48

49
```python { .api }
50
def arma2psd(A=None, B=None, rho=1., T=1., NFFT=4096, sides='default', norm=False):
51
    """
52
    Compute power spectral density given ARMA parameters.
53
    
54
    Parameters:
55
    - A: array-like, AR polynomial coefficients (None for pure MA)
56
    - B: array-like, MA polynomial coefficients (None for pure AR) 
57
    - rho: float, white noise variance
58
    - T: float, sampling period
59
    - NFFT: int, number of FFT points
60
    - sides: str, frequency range ('default', 'onesided', 'twosided')
61
    - norm: bool, normalize PSD
62
    
63
    Returns:
64
    array: Power spectral density values
65
    """
66
```
67

68
### ARMA PSD Estimation Classes
69

70
Classes for PSD estimation using ARMA and MA parameter estimation.
71

72
```python { .api }
73
class parma(data, P, Q, lag, NFFT=None, sampling=1., scale_by_freq=False):
74
    """
75
    PSD estimation using ARMA estimator.
76
    
77
    Parameters:
78
    - data: array-like, input time series
79
    - P: int, AR model order
80
    - Q: int, MA model order
81
    - lag: int, maximum lag for parameter estimation
82
    - NFFT: int, FFT length for PSD computation
83
    - sampling: float, sampling frequency
84
    - scale_by_freq: bool, scale PSD by frequency
85
    """
86

87
class pma(data, Q, M, NFFT=None, sampling=1., scale_by_freq=False):
88
    """
89
    PSD estimation using MA estimator.
90
    
91
    Parameters:
92
    - data: array-like, input time series
93
    - Q: int, MA model order
94
    - M: int, intermediate AR order for MA estimation
95
    - NFFT: int, FFT length for PSD computation
96
    - sampling: float, sampling frequency
97
    - scale_by_freq: bool, scale PSD by frequency
98
    """
99
```
100

101
## Usage Examples
102

103
### ARMA Model Parameter Estimation
104

105
```python
106
import spectrum
107
import numpy as np
108

109
# Generate ARMA(2,1) test signal
110
N = 512
111
noise = np.random.randn(N)
112

113
# True ARMA parameters
114
ar_true = [1.0, -0.5, 0.2]  # AR polynomial [1, a1, a2]
115
ma_true = [1.0, 0.3]        # MA polynomial [1, b1]
116

117
# Generate ARMA signal using scipy.signal.lfilter
118
from scipy import signal
119
arma_signal = signal.lfilter(ma_true, ar_true, noise)
120

121
# Estimate ARMA parameters
122
P, Q, lag = 2, 1, 20
123
ar_est, ma_est, var_est = spectrum.arma_estimate(arma_signal, P, Q, lag)
124

125
print(f"True AR coefficients: {ar_true[1:]}")
126
print(f"Estimated AR coefficients: {ar_est}")
127
print(f"True MA coefficients: {ma_true[1:]}")  
128
print(f"Estimated MA coefficients: {ma_est}")
129
print(f"Noise variance estimate: {var_est}")
130
```
131

132
### Pure MA Model Estimation
133

134
```python
135
import spectrum
136
import numpy as np
137

138
# Generate MA(3) signal
139
N = 256
140
noise = np.random.randn(N)
141
ma_coeffs = [1.0, 0.5, -0.3, 0.2]  # MA polynomial
142

143
from scipy import signal
144
ma_signal = signal.lfilter(ma_coeffs, [1.0], noise)
145

146
# Estimate MA parameters using long AR approximation
147
Q = 3  # True MA order
148
M = 20  # Long AR order (should be much larger than Q)
149
ma_est, var_est = spectrum.ma(ma_signal, Q, M)
150

151
print(f"True MA coefficients: {ma_coeffs[1:]}")
152
print(f"Estimated MA coefficients: {ma_est}")
153
```
154

155
### ARMA PSD Estimation and Comparison
156

157
```python
158
import spectrum
159
import numpy as np
160
import matplotlib.pyplot as plt
161

162
# Generate ARMA signal
163
N = 512
164
t = np.arange(N)
165
noise = 0.1 * np.random.randn(N)
166

167
# Create ARMA signal with known spectral characteristics
168
ar_coeffs = [1.0, -1.6, 0.64]  # Creates spectral peak
169
ma_coeffs = [1.0, 0.5]
170
from scipy import signal
171
arma_data = signal.lfilter(ma_coeffs, ar_coeffs, noise)
172

173
# ARMA PSD estimation
174
p_arma = spectrum.parma(arma_data, P=2, Q=1, lag=20)
175

176
# Compare with pure AR method
177
p_ar = spectrum.pburg(arma_data, order=10)
178

179
# Compare with non-parametric method
180
p_period = spectrum.Periodogram(arma_data)
181

182
# Plot comparison
183
plt.figure(figsize=(12, 8))
184

185
plt.subplot(2, 1, 1)
186
plt.plot(arma_data[:100])
187
plt.title('ARMA Signal (first 100 samples)')
188
plt.xlabel('Sample')
189
plt.ylabel('Amplitude')
190
plt.grid(True)
191

192
plt.subplot(2, 1, 2)
193
freqs = np.linspace(0, 0.5, len(p_arma.psd))
194
plt.semilogy(freqs, p_arma.psd, 'r-', label='ARMA(2,1)', linewidth=2)
195
plt.semilogy(freqs, p_ar.psd, 'b--', label='AR(10)', linewidth=2)
196
plt.semilogy(freqs, p_period.psd, 'g:', label='Periodogram', linewidth=2)
197
plt.xlabel('Normalized Frequency')
198
plt.ylabel('PSD')
199
plt.title('PSD Comparison: ARMA vs AR vs Periodogram')
200
plt.legend()
201
plt.grid(True)
202

203
plt.tight_layout()
204
plt.show()
205
```
206

207
### Direct ARMA to PSD Conversion
208

209
```python
210
import spectrum
211
import numpy as np
212
import matplotlib.pyplot as plt
213

214
# Define ARMA model parameters
215
ar_coeffs = np.array([1.0, -0.8, 0.15])  # AR part
216
ma_coeffs = np.array([1.0, 0.5])         # MA part
217
noise_var = 0.1
218
sampling_freq = 1.0
219

220
# Compute PSD directly from ARMA parameters
221
psd_arma = spectrum.arma2psd(A=ar_coeffs, B=ma_coeffs, rho=noise_var, 
222
                            T=1/sampling_freq, NFFT=512, sides='onesided')
223

224
# Compute pure AR PSD for comparison
225
psd_ar = spectrum.arma2psd(A=ar_coeffs, B=None, rho=noise_var, 
226
                          T=1/sampling_freq, NFFT=512, sides='onesided')
227

228
# Compute pure MA PSD for comparison  
229
psd_ma = spectrum.arma2psd(A=None, B=ma_coeffs, rho=noise_var, 
230
                          T=1/sampling_freq, NFFT=512, sides='onesided')
231

232
# Plot theoretical PSDs
233
freqs = np.linspace(0, sampling_freq/2, len(psd_arma))
234
plt.figure(figsize=(10, 6))
235
plt.semilogy(freqs, psd_arma, 'r-', label='ARMA(2,1)', linewidth=2)
236
plt.semilogy(freqs, psd_ar, 'b--', label='AR(2)', linewidth=2)
237
plt.semilogy(freqs, psd_ma, 'g:', label='MA(1)', linewidth=2)
238
plt.xlabel('Frequency (Hz)')
239
plt.ylabel('PSD')
240
plt.title('Theoretical ARMA PSD Components')
241
plt.legend()
242
plt.grid(True)
243
plt.show()
244
```