Sites about Auto Regressive Moving Average (ARMA)

Since I need to design a controller for the Adaptive optics, there is a topic about Kalman filtering and Optimal Control. Here is a small collection of helpful materials I have found so far. Introduction to Kalman filters is here. For an excellent web site, see Welch/Bishop's KF page. Some tutorials, references, and research related to the Kalman filter.Kalman filter Toolbox can be reached here.




MATLAB functions for ARMA
Maximum Likelihood Estimation
This section explains how the garchfit estimation engine uses maximum likelihood to estimate the parameters needed to fit the specified models to a given univariate return series

Lots of MATLAB m-files for ARMA are in:

• homeworks pages.
• ARMA for control
• pattern recognition
• Kalman filtering
• Oxford works in Kalman filters
  

Discussions

MATLAB code for fitting ARMA/ARIMA models?

Does anyone out there have a MATLAB code for fitting ARMA models (with
specified autoregressive order p and moving average order q) to time
series data?

>> help stats

Statistics Toolbox.
Version 3.0 (R12) 1-Sep-2000

see System Identification TB

Please see "ARMAX", and "PEM" in the System Identification toolbox.

For time series, the function "PRONY" in Signal Processing toolbox, and the
function "HPRONY" (which is based on PRONY) in the HOSA (Higher Order
Spectral Analysis) toolbox can also be used.



u=iddata(timeseries)
m = armax(u,[p q]) %ARMA(p,q)

result:
Discrete-time IDPOLY model: A(q)y(t) = C(q)e(t)
A(q) = 1 - 1.216 q^-1 + 0.7781 q^-2

C(q) = 1 - 1.362 q^-1 + 0.8845 q^-2 + 0.09506 q^-3

Estimated using ARMAX from data set z
Loss function 0.000768896 and FPE 0.00084009
Sampling interval: 1

Is there anyone have the idea how to use matlab code to fitting it?

Thanks.


ARIMA model can be created by using differenced data with ARMAX. For
example, use diff(x) rather than x as output data and then using ARMAX
command. A pole at 1 can be added to resulting model to achieve a true ARIMA
model.

m = armax(diff(x),[p,q])
m.a = conv(m.a,[1 -1]);


Also note that you can specify a differencing order using DIFF(X,ORDER) based on the trend (stationarity) observed in your data.

And remember too if you do forecasting this way that you have to integrate the resulting predicted time series by recursively summing the predicted values using the last observation value as your initial condition.

ARMA estimation [repost]

MATLAB has built in functions in its toolbox to estimate ARMA/GARCH
parameteres once a model is specified (garchfit). It also tells me that it is
using Maximum Likelihood Estimataion procedure, and it seems to do it
very fast. I am curious, if anybody knows such detail, whether it is
doing an MLE estimation with i.i.d. assumption (which is clearly
false), or does it actually do the recursive factorization of joint
pdf?

I think you should check the info pages on the Mathworks
website about the Econometrics Toolbox. You can download
pdf's with a lot of information about the toolbox and the
computational details. I was curious...it seems to be
a nonlinear programming approach with iterative
nature.

This section of the documentation explains the MLE calculation:

http://www.mathworks.com/help/toolbox/econ/f1-80052.html

ARMA/ARIMA model identification

I am trying to use the System Identification Toolbox to determin the best ARMA or ARIMA model to use. I am getting errors from the arxstruc and selctruc. I also get errors if I specify the ARMA model with out a 1 for the MA order.

There are errors in the docs which say conflicting things about which models forms will work (e.g. "Note Use selstruc for single-output systems only. selstruc supports both single-input and multiple-input systems.")

AR(I)MA sytem identification is inherently difficult.
Don't expect anything such routine to work unconditionally.

ARMA model fitting

I've been working in developing input/output ARMA model for surface roughness profile. I am using MATLAB armax command for that, which has a flexibility to choose my model's order. I was just wondering, is there any way I can plot the residuals, obtained from the ARMA model. My intension is to verify the adequacy of the developed model, that's why thinking to figure out the aucorrelation function of residuals. But I'm really confused regarding how to do that.

obtain the numerator and denominator of the ARMA model (should be in the ARMA model structure), then use filter() with these parameters along with the ARMA system input as inputs. Take the output from filter() and compare it to the original output to see residuals and all that.

Thank you very much for your suggestion. The thing is MATLAB returns ARMA model as an IDPOLY output. For an instance, for my case the intended structure is like,
Discrete-time IDPOLY model:
A(q)y(t) = B(q)u(t) + C(q)e(t)
where, y(t)=> output, u(t)=> input
And
A(q) = 1 - 0.8356 q^-1 + 0.1549 q^-2

B(q) = 0.08899 q^-2 + 0.08928 q^-3

C(q) = 1 + 0.2435 q^-1 + 0.3926 q^-2

For this kind of model, what should be the numerator and denominator? It would be great, if you state a bit more explicitly.



Hi Hafiz,
As a difference equation your model is:

y(t)=0.8356*y(t-1)-0.1549*y(t-2)+0.08899*u(t-2)+0.08928*u(t-3)+e(t)+0.2435*e(t-1)+0.3926*e(t-2)

y(t) is the output, u(t) is the exogenous term, and e(t) is the white noise disturbance. In this model, you are writing the output as a sum of two convolutions, the first one being a convolution of the input, u(t), with the impulse response, and the 2nd convolution being a filtering of the disturbance, or noise term, e(t). Hence, the transfer function is the ratio of:

(0.08899*z^{-2}+0.08928*z^{-3})/(1-0.8356*z^{-1}+0.1549*z^{-2})

To view:

fvtool([0 0 0.08899 0.08928],[1 -0.8356 0.1549])

Hope that helps.


Look up functions RESID, COMPARE, PREDICT and PE (in System Identification
Toolbox). In particular, RESID will plot the autocorrelation of residuals.


I need something to know about 'resid' command. I'm using it like this:

d= resid (arma_model, DAT, 'corr');
plot (d);

The plots return the autocorrelation of error and cross correlation between the error and the input. I'm wondering is there any way I can plot the autocorrelation of the residuals in normalized form (in a scale of y between 0 to 1)?

arma/armax

I'm having some issues with the armax function. I have an single input-output time series dataset, and all I want to do is extimate the parameters of an arma(2,3) model, without the exogenous term. How do I go about this?

Hi Brian,

You can do

% creating some model data ARMA(1,2)
A=[1 -0.9];
B = [ 1 0.2000 0.4000];
reset(RandStream.getDefaultStream);
y = filter(B,A,randn(1000,1));
Data = iddata(y,[],1); % no exogenous term
Model =armax(Data,[1 2]);
% displaying Model

Model
Discrete-time IDPOLY model: A(q)y(t) = C(q)e(t)
A(q) = 1 - 0.9012 q^-1

C(q) = 1 + 0.2259 q^-1 + 0.4615 q^-2

Estimated using ARMAX on data set Data
Loss function 0.993907 and FPE 0.999895
Sampling interval: 1

Compare A(q) to A and C(q) to B

Hope that helps.







MATLAB software about ARMA modelling

ARMA PARAMETER ESTIMATION

This function provides an ARMA spectral estimate which is maximum entropy satisfying correlation constraint (number of poles)and cepstrum
constrains (number of ceros)

The function requires 3 inputs: Input signal, Order of denonominator, Order of Numerator and the output variable are: numerator coeffients, denominator coeficients and square-root of input noise power.

by Miguel Angel Lagunas Hernandez

Estimating the parameters of AR-MA sequences is fundamental. Here we present another method for ARM. Download Now

Automatic Spectral Analysis

by Stijn de Waele

03 Jul 2003 (Updated 09 May 2005)

No BSD License

Automatic spectral analysis for irregular sampling/missing data, analysis of spectral subband.

Download Now

Accurate estimates of the autocorrelation or power spectrum can be obtained with a parametric model (AR, MA or ARMA). With automatic inference, not only the model parameters but also the model structure are determined from the data. It is assumed that the ARMASA toolbox is present. This toolbox can be downloaded from the MATLAB Central file exchange at www.mathworks.com
The applications of this toolbox are:
- Reduced statistics ARMAsel: A compact yet accurate ARMA model is obtained based on a given power spectrum. Can be used for generation of colored noise with a prescribed spectrum.
- ARfil algorithm: The analysis of missing data/irregularly sampled signals
- Subband analysis: Accurate analysis of a part of the power spectrum
- Vector Autoregressive modeling: The automatic analysis of auto- and crosscorrelations and spectra
- Detection: Generally applicable test statistic to determine whether two signals have been generated by the same process or not. Based on the Kullback-Leibler index or Likelihood Ratio.
- Analysis of segments of data, possibly of unequal length.

For background information see my PhD thesis, available at http://www.dcsc.tudelft.nl/Research/PubSSC/thesis_sdewaele.html

MATLAB parallelisation for the Adaptive Optics Simulators

Recently in the article[1] there were proposed interesting ideas about parallelisation of MATLAB. In Euro50 project of the Extremely Large Telescope, the model of the entire telescope is built on the MATLAB and calculations are obviously beyond the capabilities of personal computers.

MATLAB does not include any parallel functionality however there are basic interfaces to the TCP/IP network stack, C and Fortran. By building upon these interfaces there have been of the order of 30 attempts (as in the survey R. Choy, Parallel MATLAB survey) to produce ``toolkits'' to allow MATLAB to be used in a parallel fashion. The approaches of communications between nodes that were used in those toolkits either take a low-level approach using commands similar to Message Passing Interface (MPI) libraries or use a more high level approach choosing to use simpler commands that resemble more closely those of MATLAB. The problem with MPI is that it is difficult to use with MATLAB[1].

According to[1], the acceleration of the simulation's execution using compilation to MEX is significant: experiment showed that MATLAB could indeed compete with the traditional HPC languages in raw performance, typically attaining 90% of the performance of Fortran 90 on calculations.

As a starting point, the parallelisation toolkit for MATLAB written by Einar Heiberg was used.
The Matlab parallelization toolkit is released as Open Source, uses a Master/Slave paradigm and is most suitable to problems where the amount of communication is low. The toolkit by D.Moraru called MatlabWS that is used for Euro50 is not available publicly. But the results of the investigations are interesting. The authors of[1] found out that in practise, a bigger problem is not a network bandwidth but a latency introduced by network card driver. As an example, using the Heiberg toolkit have shown that on 100Mbps Ethernet there is minimum period of 35ms of latency involved in any communications between MATLAB instances on separate cluster nodes regardless of message size. While moving to gigabit Ethernet would increase bandwidth, it would have little impact on latency[1].

Another approach is MPITB written by Javier Fernandez Baldomero that uses MPI. According to the webpage, PC MATLAB Linux users in a cluster with several PCs can use MPITB in order to call MPI library routines from within the MATLAB interpreter. Depending on your licensing scheme (node/user-based), additional licenses might be required to spawn additional MATLAB processes on other cluster nodes. Currently processes can be spawned and arranged in topologies, MATLAB variables can be sent/received.

As a conclusion, the improvements in latency have resulted in a reduction in typical model run time from 70 hours to 24 hours. It is hoped that architectural changes that better exploit a lower latency environment can reduce this still further. In addition another toolkit MPITB offers comparable latencies over Ethernet, using native MPI. MPITB uses the LAM MPI implementation. However this toolkit is more complex to configure and use. MatlabWS's combination of the ease of use of the Heiberg toolkit with the performance of MPITB make it a very compelling product.

[TBD] Ideas about parallelisation of GNU/Octave are collected here.

References:
[1] Browne, M., Andersen, T., Enmark,
A., Moraru, D., and Shearer, A., "Parallelization of MATLAB for Euro50 Integrated
Modeling", Proc. SPIE, Vol. 5497, 2004.

MATLAB vs Octave Benchmark

The benchmark is performed by a script that is a part of AORTA simulator. For the time-consumption measure, standard tic-toc combination is used. For the memory consumption, UNIX command is executed during the computation to obtain actual memory consumption.

The hardware used was Dell Latitude E5100 laptop with Intel Celeron 2.2GHz processor, 2Gb DDR RAM, and ATA-100 160Gb HDD. The test was run in a system's terminal (graphical X Window system was shut down) so the interference with other software was negligible.

Since the simulator uses Fast Fourier Transform (FFT) and matrix operations heavily, the benchmark consists of those operations. In order to gain maximum performance, all matrices were power of two: from 2x2 up to 8192x8192 .

Those numbers, of course, are not absolute and change from version to version. But the according to our results, the tendency remains the same: Octave tends to be faster on Fourier transform than MATLAB but consumes more memory.

Comparison of computation time

The basic matrix matrix operations, such as Gauss and point-wise multiplication as well as summation, the performance of MATLAB and Octave are approximately the same. The differences are in computational time of Fourier transform. As one can see from the Fig 1.1, Octave outperforms MATLAB on FFT transforms on any matrix sizes.

a) b) Figure 1.1: Median computation time versus matrix size for: a) FFT function and for b) FFTSHIFT function.

For example, on typical size of matrices in adaptive optics simulations (thousands by thousands of elements), Octave is almost 3 times faster than MATLAB with similar memory consumption. However for 8192x8192 matrix, the ``memory exhausted or requested size too large for range of Octave's index type'' message appeared and there is no datapoint for such a large matrix. It is also notable that on Windows and MATLAB v2009 an attempt to calculate huge matrices (more than 5000-6000 elements in each direction) results to the same ``out of memory'' error message.

Hence, except matrices that are bigger than 8192x8192, GNU/Octave is almost 3 times faster than MATLAB in FFT.

Comparison of memory consumption

Memory consumption for Octave is generally larger than for MATLAB as it is apparent from the Fig. 1.2. For small matrices (less than 2048x2048) it can be 2.5-3 times more memory, and for large matrices ( and more) that is typically 10-20%. Consequently, on the available hardware with 2Gb RAM the calculation of FFT of matrix 8192x8192 elements is not possible.

Figure 1.2: Median memory consumption versus matrix size for: a) FFT function and for b) FFTSHIFT function.

a) b)

Although Octave consumes considerably more memory than MATLAB while performing FFT and thus may run out of memory on huge matrices.

Routh array, Routh Criterion and how to calculate the Routh Hurwitz array

The post is available also in PDF format here[mirror].

Trying to find a reasonable explanation of how to calculate Routh array, I found one here (while in book are usually given cryptic and tricky ways, even in Goodwin-Graebe).

Routh criterion

The Routh-Hurwitz stability criterion is a necessary (and frequently sufficient) method to establish the stability of a single-input, single-output (SISO), linear time invariant (LTI) control system.

One of the most popular algorithms to determine whether or not a polynomial is strictly Hurwitz, is Routh's algorithm. Consider a polynomial and its associated Routh array (see below). Then the number of roots with real part greater than zero is equal to the number of sign changes in the first column of the array.

Routh array

A tabular method of Routh criterion can be used to determine the stability when the roots of a higher order characteristic polynomial are difficult to obtain. Consider the polynomial . The first two rows (two highest ones) of the Routh's array are formed using coefficients in the increasing order:

Row	Power	First column	Second column	Third column
n				...
n-1				...

If there are no more coefficients then just add zeros in the table.

Then for each row after third row we need to calculate , , and so on. Now starting from row 3, each row contains the b's, c's, d's and so on as follows:

Row	Power	First column	Second column	Third column
n				...
n-1				...
n-2				...
n-3				...
n-4				...
&vellip#vdots;
0			0	0

Now we need to calculate the coefficients and so on. That can be done using the following determinants formulas:

and so on until . Remember to add zeros in the table when no more coefficients left - it will alleviate the calculations and make them more straightforward. The principle is the following: when one calculates , the first column in the determinant is always [column_1; row_1, row_2]"/> of the Routh table, e.g. .

Second column in the determinant contains , e.g. for calculations of .

The following illustration shows the principle:

Here dark-green rectangle depicts the first column of the determinant in the calculations of , dark-green circle - denominator in the equation for . Bright-green rectangle depicts the second column in the determinant for calculations of . In the same fashion, dark-blue rectangle depicts the first column of the determinant in the calculations of , dark-blue circle - denominator in the equation for . Bright-blue rectangle depicts the second column in the determinant for calculations of , and so on.

Example of calculation of Routh array

Consider the polynomial . As above, the first two rows of the Routh's array are formed using coefficients in the increasing order:

Row	Power	First column	Second column	Third column
n		1	3	1
n-1		1	2	0 ...

Now we need to calculate using mentioned above equations.

Thus our table will look like:

Row	Power	First column	Second column	Third column
n		1	3	1
n-1		1	2	0 ...
n-2		1	1	0 ...
n-3		1	0	0 ...
n-4		1	0	0 ...

From the array we note that there are no sign changes in the first column. According to Routh's criterion, this means that is a strictly Hurwitz polynomial.