Comparing SISO ARX one step ahead prediction to inf step ahead prediction (simulation). This uses this branch of DTK: https://github.com/moorepants/DynamicistToolKit/tree/system-id

prediction_simulation.py

import numpy as np
from scipy.optimize import leastsq
import pandas
import matplotlib.pyplot as plt
from dtk.process import butterworth
from dtk.systemid import ARX

df = pandas.read_csv('RandomSlowSpeed1.txt', delimiter='\t')
sample_rate = 100.0  # hz
cutoff = 6.0  # hz
df['Filtered Belt Speed'] = butterworth(df['Belt Velocity (m/s)'].values,
                                        cutoff, sample_rate)
df['Filtered Pitch Moment'] = butterworth(df['Pitch Moment (Nm)'].values,
                                          cutoff, sample_rate)
df['Belt Acceleration'] = (df['Filtered Belt Speed'].diff() /
                           df['Timestamp (s)'].diff())
df['Belt Acceleration'] = df['Belt Acceleration'].fillna(0.0)

# subtract means
df['Belt Acceleration'] = df['Belt Acceleration'] - df['Belt Acceleration'].mean()
df['Filtered Pitch Moment'] = df['Filtered Pitch Moment'] - df['Filtered Pitch Moment'].mean()

arx = ARX(df['Belt Acceleration'].values,
          df['Filtered Pitch Moment'].values,
          na=2, nb=3, nk=0)

arx.form_regressor()
arx.form_a_b()
arx.solve()
N = 1000
x0 = arx.find_initial_condition(df['Belt Acceleration'][:N].values,
                                guess=np.array([-1.52968045,  1.18207912,  1.30084501]))
y_sim = arx.simulate(df['Belt Acceleration'][:N].values, x0)

# Now find the answer through brute force optimization with discrete
# simulation.

def simulate(theta, u):
    y = np.zeros_like(u)
    for i in range(len(y)):
        if i >= 2:
            y[i] = (theta[0] * -y[i - 2] +
                    theta[1] * -y[i - 1] +
                    theta[2] * u[i - 2] +
                    theta[3] * u[i - 1] +
                    theta[4] * u[i])
    return y

def residuals(theta, u, y):
    return simulate(theta, u) - y

result = leastsq(residuals, arx.solution, args=(df['Belt Acceleration'][:N].values,
                                                df['Filtered Pitch Moment'][:N].values))

y_sim2 = simulate(result[0], df['Belt Acceleration'][:N].values)

arx.solution = result[0]
y_sim3 = arx.simulate(df['Belt Acceleration'][:N].values, arx.y[:3])

plt.plot(df['Timestamp (s)'][:len(y_sim)], arx.y[:len(y_sim)], 'k-')

ax = plt.gca()
y_lim = ax.get_ylim()

plt.plot(df['Timestamp (s)'][:len(y_sim)], y_sim, 'b')
plt.plot(df['Timestamp (s)'][:len(y_sim)], y_sim2, 'r')
plt.plot(df['Timestamp (s)'][:len(y_sim)], y_sim3, 'g')

plt.ylim(y_lim)
#plt.legend(('Measured', 'Simulated {:1.1%}'.format(r2)))

plt.show()