ARMA Methods

def arma_estimate(X, P, Q, lag):
    """
    Estimate ARMA(P,Q) parameters using modified Yule-Walker technique.
    
    Parameters:
    - X: array-like, input time series data
    - P: int, AR model order
    - Q: int, MA model order  
    - lag: int, maximum lag for correlation computation
    
    Returns:
    tuple: (AR parameters array, MA parameters array, white noise variance float)
    """

def ma(X, Q, M):
    """
    Moving average estimator using long AR model approach.
    
    Parameters:
    - X: array-like, input time series data
    - Q: int, MA model order
    - M: int, intermediate AR model order (should be >> Q)
    
    Returns:
    tuple: (MA parameters array, white noise variance estimate float)
    """

def arma2psd(A=None, B=None, rho=1., T=1., NFFT=4096, sides='default', norm=False):
    """
    Compute power spectral density given ARMA parameters.
    
    Parameters:
    - A: array-like, AR polynomial coefficients (None for pure MA)
    - B: array-like, MA polynomial coefficients (None for pure AR) 
    - rho: float, white noise variance
    - T: float, sampling period
    - NFFT: int, number of FFT points
    - sides: str, frequency range ('default', 'onesided', 'twosided')
    - norm: bool, normalize PSD
    
    Returns:
    array: Power spectral density values
    """

class parma(data, P, Q, lag, NFFT=None, sampling=1., scale_by_freq=False):
    """
    PSD estimation using ARMA estimator.
    
    Parameters:
    - data: array-like, input time series
    - P: int, AR model order
    - Q: int, MA model order
    - lag: int, maximum lag for parameter estimation
    - NFFT: int, FFT length for PSD computation
    - sampling: float, sampling frequency
    - scale_by_freq: bool, scale PSD by frequency
    """

class pma(data, Q, M, NFFT=None, sampling=1., scale_by_freq=False):
    """
    PSD estimation using MA estimator.
    
    Parameters:
    - data: array-like, input time series
    - Q: int, MA model order
    - M: int, intermediate AR order for MA estimation
    - NFFT: int, FFT length for PSD computation
    - sampling: float, sampling frequency
    - scale_by_freq: bool, scale PSD by frequency
    """

import spectrum
import numpy as np

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

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

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

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

print(f"True AR coefficients: {ar_true[1:]}")
print(f"Estimated AR coefficients: {ar_est}")
print(f"True MA coefficients: {ma_true[1:]}")  
print(f"Estimated MA coefficients: {ma_est}")
print(f"Noise variance estimate: {var_est}")

import spectrum
import numpy as np

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

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

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

print(f"True MA coefficients: {ma_coeffs[1:]}")
print(f"Estimated MA coefficients: {ma_est}")

import spectrum
import numpy as np
import matplotlib.pyplot as plt

# Generate ARMA signal
N = 512
t = np.arange(N)
noise = 0.1 * np.random.randn(N)

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

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

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

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

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

plt.subplot(2, 1, 1)
plt.plot(arma_data[:100])
plt.title('ARMA Signal (first 100 samples)')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)

plt.subplot(2, 1, 2)
freqs = np.linspace(0, 0.5, len(p_arma.psd))
plt.semilogy(freqs, p_arma.psd, 'r-', label='ARMA(2,1)', linewidth=2)
plt.semilogy(freqs, p_ar.psd, 'b--', label='AR(10)', linewidth=2)
plt.semilogy(freqs, p_period.psd, 'g:', label='Periodogram', linewidth=2)
plt.xlabel('Normalized Frequency')
plt.ylabel('PSD')
plt.title('PSD Comparison: ARMA vs AR vs Periodogram')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

import spectrum
import numpy as np
import matplotlib.pyplot as plt

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

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

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

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

# Plot theoretical PSDs
freqs = np.linspace(0, sampling_freq/2, len(psd_arma))
plt.figure(figsize=(10, 6))
plt.semilogy(freqs, psd_arma, 'r-', label='ARMA(2,1)', linewidth=2)
plt.semilogy(freqs, psd_ar, 'b--', label='AR(2)', linewidth=2)
plt.semilogy(freqs, psd_ma, 'g:', label='MA(1)', linewidth=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.title('Theoretical ARMA PSD Components')
plt.legend()
plt.grid(True)
plt.show()