#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# Doc/
# Doc/CHANGES
# Doc/README
# Doc/index.css
# Doc/index.html
# Matlab/
# Matlab/Sp/
# Matlab/Sp/Drivers/
# Matlab/Sp/Drivers/ardem.m
# Matlab/Sp/Src/
# Matlab/Sp/Src/acf.m
# Matlab/Sp/Src/adjph.m
# Matlab/Sp/Src/arconf.m
# Matlab/Sp/Src/arfit.m
# Matlab/Sp/Src/armode.m
# Matlab/Sp/Src/arord.m
# Matlab/Sp/Src/arqr.m
# Matlab/Sp/Src/arres.m
# Matlab/Sp/Src/arsim.m
# Matlab/Sp/Src/tquant.m
# This archive created: Wed Sep 26 11:59:36 2001
export PATH; PATH=/bin:$PATH
if test ! -d 'Doc'
then
mkdir 'Doc'
fi
cd 'Doc'
if test -f 'CHANGES'
then
echo shar: will not over-write existing file "'CHANGES'"
else
cat << "SHAR_EOF" > 'CHANGES'
14-Oct-00: * Minor revision of the papers describing the algorithms.
* Minor changes to some modules to improve error handling and
warning messages.
06-Sep-00: * The ARfit homepage has moved to www.math.nyu.edu/~tapio/arfit/.
15-Apr-00: * Changed calling sequence of Matlab function FZERO in
TQUANT to old format in order to ensure compatibility
with Windows/Mac versions of Matlab. [Suggested by
Vico Klump]
10-Jan-00: * The two papers describing the algorithms implemented in ARfit
have been extensively revised.
* Confidence intervals are now based on Student's t
distribution, no longer on the chi-squared distribution.
95% margins of error are returned instead of standard
errors.
* The regularization in ARQR has changed.
* The margin of error for the period of a purely relaxatory mode
with infinite period is now zero, not NaN as before.
17-Dec-99: * Renamed the demo-module from ARFIT to ARDEM.
* Renamed the least squares estimation module from AR to
ARFIT (in order to avoid conflicts with ar.m in the system
identification toolbox).
* Rewrote some of the help entries and annotations.
* Changed output of ARCONF. The approximate confidence intervals
are now returned as separate variables, no longer as imaginary
parts of the parameter matrices.
* The modification MSC of Schwarz's Bayesian Criterion for order
selection is no longer supported.
* Changed definition of the excitation (change affects only
higher-order models).
12-Aug-98: * Fixed bug (in AR) that affected the computation of confidence
intervals for the intercept vector. (The condition-improving
scaling was not undone before the matrix Uinv was returned.)
02-Aug-98: * Fixed bug in ARSIM. (ARSIM did not simulate multivariate AR(1)
processes correctly.)
* Replaced DOT(a,b) by SUM(a.*b) because DOT(a,b) is neither
part of older Matlab distributions nor part of Octave.
Only ADJPH was affected by this change, which increases the
portability of ARfit.
* Forced dot_lam in ARMODE to be a column vector in order to
ensure compatibility with Octave.
* Removed (from AR) calls of the function LOWER in order to
ensure compatibility with Octave.
* Readability of the code has been improved.
* All modules have been tested under Matlab 5.
06-Sep-97: * Release of ARfit version 2.0 with accelerated
computation of order selection criteria.
SHAR_EOF
fi # end of overwriting check
if test -f 'README'
then
echo shar: will not over-write existing file "'README'"
else
cat << "SHAR_EOF" > 'README'
%% arfit.tar
%% Authors: Tapio Schneider and Arnold Neumaier
%% Contents:
CHANGES A history of recent revisions of the programs
schneider1.ps PostScript of "Estimation of parameters and eigenmodes of
multivariate autoregressive models"
schneider1.tex LaTeX source of schneider1.ps (in ACM house style)
schneider2.ps PostScript of "Algorithm XXX: ARfit - A Matlab package for
the estimation of parameters and eigenmodes of multivariate
autoregressive models"
schneider2.tex LaTeX source of schneider2.ps (in ACM house style)
index.html HTML page with a short description of how to install and use
index.css the Matlab modules.
acf.m |
adjph.m |
arconf.m |
ardem.m |
arfit.m | Matlab modules of the ARfit package (see index.html for
armode.m | descriptions of the modules)
arord.m |
arqr.m |
arres.m |
arsim.m |
tquant.m |
SHAR_EOF
fi # end of overwriting check
if test -f 'index.css'
then
echo shar: will not over-write existing file "'index.css'"
else
cat << "SHAR_EOF" > 'index.css'
/*
* Style sheet for software pages
*/
body{
font-family: Times;
font-size: medium;
background-color: white;
margin-left: 4em;
width: 35em
}
td{
font-family: Times;
font-size: medium
}
p.ref{ /* for bibliography */
margin-left: 4em;
text-indent: -2.5em
}
code{
font-family: Courier, monospace
}
p.code{
margin-left: 2.5em;
text-align: left;
font-family: Courier, monospace
}
h1{
font-size: large;
text-align: center;
}
h2{
text-align: left;
margin-top: 1.5ex
}
div.cent{
text-align: center
}
SHAR_EOF
fi # end of overwriting check
if test -f 'index.html'
then
echo shar: will not over-write existing file "'index.html'"
else
cat << "SHAR_EOF" > 'index.html'
ARfit: Fitting Multivariate Autoregressive Models
ARfit: A Matlab package for the estimation of parameters and
eigenmodes of multivariate autoregressive models
estimating parameters of multivariate autoregressive (AR)
models,
diagnostic checking of fitted AR models, and
analyzing eigenmodes of fitted AR models.
The algorithms implemented in ARfit are described in the following
papers:
A. Neumaier and T. Schneider, 2000: Estimation of
parameters and eigenmodes of multivariate autoregressive
models. To appear in ACM Trans. Math. Softw.
T. Schneider and A. Neumaier, 2000: Algorithm:
ARfit - A Matlab package for the estimation of parameters and
eigenmodes of multivariate autoregressive models. To appear in
ACM Trans. Math. Softw.
ARfit has been successfully tested under Matlab 3 and later
versions, including Matlab 5.
The ARfit package consists of a number of Matlab modules and the
file CHANGES with a history of recent
revisions of the programs.
To install ARfit, copy the Matlab modules (files with names
ending in .m) into a directory that is accessible by
Matlab. Starting Matlab and invoking Matlab's online help
function
help filename
calls up detailed information on the purpose and the calling
syntax of the module filename.m. The script
ardem.m demonstrates the basic features of the modules contained
in ARfit.
Multiplies a complex vector by a phase factor such that the
real part and the imaginary part of the vector are orthogonal
and the norm of the real part is greater than or equal to the
norm of the imaginary part. ADJPH is required by ARMODE to
normalize the eigenmodes of an AR model.
Eigendecomposition of AR model. For a fitted AR model,
ARMODE computes eigenmodes and their associated oscillation
periods and damping times, as well as approximate confidence
intervals for the eigenmodes, periods, and damping times.
Diagnostic checking of the residuals of a fitted
model. Computes the time series of residuals. The modified
multivariate portmanteau statistic of Li & McLeod (1981) is
used to test the residuals for uncorrelatedness.
Tapio
Schneider
Courant Institute of Mathematical Sciences
New York University
251 Mercer Street
New York, NY 10012
Arnold
Neumaier
Institut für Mathematik
Universität Wien
A-1090 Wien
Austria
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Matlab'
then
mkdir 'Matlab'
fi
cd 'Matlab'
if test ! -d 'Sp'
then
mkdir 'Sp'
fi
cd 'Sp'
if test ! -d 'Drivers'
then
mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'ardem.m'
then
echo shar: will not over-write existing file "'ardem.m'"
else
cat << "SHAR_EOF" > 'ardem.m'
%ARDEM Demonstrates modules of the ARfit package.
% Revised: 30-Dec-99 Tapio Schneider
format short
format compact
echo on
clc
% ARfit is a collection of Matlab modules for the modeling of
% multivariate time series with autoregressive (AR) models. ARfit
% contains modules for estimating parameters of AR models from given
% time series data; for checking the adequacy of an estimated
% AR model; for analyzing eigenmodes of an estimated AR model; and
% for simulating AR processes.
%
% This demo illustrates the use of ARfit with a bivariate AR(2)
% process
%
% v(k,:)' = w' + A1*v(k-1,:)' + A2*v(k-2,:)' + eta(k,:)',
%
% where the row vectors eta(k,:) are independent and identically
% distributed Gaussian noise vectors with zero mean and covariance
% matrix C. The kth row v(k,:) of the 2-column matrix v represents an
% observation of the process at instant k. The intercept vector w is
% included to allow for a nonzero mean of the AR(p) process.
pause % Press any key to continue.
clc
% Let us simulate observations from a bivariate AR(2) process,
% choosing the parameters
w = [ 0.25 ; 0.1 ];
% for the intercept vector,
A1 = [ 0.4 1.2; 0.3 0.7 ];
% and
A2 = [ 0.35 -0.3; -0.4 -0.5 ];
% for the AR coefficient matrices, and
C = [ 1.00 0.50; 0.50 1.50 ];
% for the noise covariance matrix. The two 2x2 matrices A1 and A2 are
% assembled into a single 2x4 coefficient matrix:
A = [ A1 A2 ];
% We use the module ARSIM to simulate 200 observations of this AR
% process:
v = arsim(w, A, C, 200);
pause % Press any key to continue.
clc
% Suppose that we have no information about how the time series
% data v are generated, but we want to try to fit an AR model to the
% time series. That is, we must estimate the AR model parameters,
% including the model order. Assuming that the correct model order
% lies between 1 and 5, we use the module ARFIT to determine the
% optimum model order using Schwarz's Bayesian Criterion (SBC):
pmin = 1;
pmax = 5;
[west, Aest, Cest, SBC, FPE, th] = arfit(v, pmin, pmax);
% The output arguments west, Aest, and Cest of ARFIT are the
% estimates of the intercept vector w, of the coefficient matrix A,
% and of the noise covariance matrix C. (The matrix th will be needed
% later in the computation of confidence intervals.) These parameters
% are estimated for a model of order popt, where the optimum model
% order popt is chosen as the minimizer of an approximation to
% Schwarz's Bayesian Criterion. The selected order popt in our
% example is:
m = 2; % dimension of the state space
popt = size(Aest, 2) / m;
disp(['popt = ', num2str(popt)])
pause % Press any key to continue.
clc
% Besides the parameter estimates for the selected model, ARFIT has
% also returned approximations to Schwarz's Bayesian Criterion SBC
% and to Akaike's Final Prediction Error FPE, each for models of
% order pmin:pmax. In this demo, the model order was chosen as the
% minimizer of SBC. Here are the SBCs for the fitted models of order
% 1,...,5:
disp(SBC)
% To see if using Akaike's Final Prediction Error as a criterion to
% select the model order would have resulted in the same optimum
% model order, compare the FPEs:
disp(FPE)
% Employing FPE as the order selection criterion, the optimum model
% order would have been chosen as the minimizer of FPE. The values of
% the order selection criteria are approximations in that in
% evaluating an order selection criterion for a model of order p <
% pmax, the first pmax-p observations of the time series are ignored.
pause % Press any key to continue.
clc
% Next it is necessary to check whether the fitted model is adequate
% to represent the time series v. A necessary condition for model
% adequacy is the uncorrelatedness of the residuals. The module ARRES
[siglev,res] = arres(west,Aest,v);
% returns the time series of residuals res as well as the
% significance level siglev with which a modified Li-McLeod
% portmanteau test rejects the null hypothesis that the residuals are
% uncorrelated. A model passes this test if, say, siglev > 0.05. In
% our example, the significance level of the modified Li-McLeod
% portmanteau statistic is
disp(siglev);
% (If siglev is persistently smaller than about 0.05, even over
% several runs of this demo, then there is most likely a problem with
% the random number generator that ARSIM used in the simulation of
% v.)
pause % Press any key to continue.
clc
if exist('xcorr') % If the Signal Processing Toolbox is installed, plot
% autocorrelation function of residuals ...
% Using ACF, one can also plot the autocorrelation function of, say,
% the first component of the time series of residuals:
acf(res(:,1));
% 95% of the autocorrelations for lag > 0 should lie between the
% dashdotted confidence limits for the autocorrelations of an IID
% process of the same length as the time series of residuals res.
% Since uncorrelatedness of the residuals is only a necessary
% condition for model adequacy, further diagnostic tests should, in
% practice, be performed. However, we shall go on to the estimation
% of confidence intervals.
pause % Press any key to continue.
end
clc
% Being reasonably confident that the model adequately represents the
% data, we compute confidence intervals for the AR parameters with
% the module ARCONF:
[Aerr, werr] = arconf(Aest, Cest, west, th);
% The output arguments Aerr and werr contain margins of error such
% that (Aest +/- Aerr) and (west +/- werr) are approximate 95%
% confidence intervals for the individual AR coefficients and for the
% components of the intercept vector w. Here is the estimated
% intercept vector with margins of error for the individual
% components in the second column:
disp([west werr])
% For comparison, the `true' vector as used in the simulation:
disp(w)
pause % Press any key to display the other parameter estimates.
clc
% The estimated coefficient matrix:
disp(Aest)
% with margins of error for its elements:
disp(Aerr)
% For comparison, the `true' AR coefficients:
disp(A)
pause % Press any key to continue.
echo off
% Compute `true' eigenmodes from model parameters:
% Eigenvectors and eigenvalues of corresponding AR(1) coefficient
% matrix:
[Strue,EvTrue] = eig([A; eye(2) zeros(2)]);
EvTrue = diag(EvTrue)'; % `true' eigenvalues
Strue = adjph(Strue); % adjust phase of eigenmodes
clc
echo on
% Finally, ARMODE computes the eigendecomposition of the fitted AR
% model:
[S, Serr, per, tau, exctn] = armode(Aest, Cest, th);
% The columns of S are the estimated eigenmodes:
disp(S)
% with margins of error Serr:
disp(Serr)
% The intervals (S +/- Serr) are approximate 95% confidence intervals
% for the individual components of the eigenmodes. Compare the
% estimated eigenmodes above with the eigenmodes obtained from the
% `true' model parameters:
disp(Strue(3:4,:))
% (Note that the estimated modes can be a permutation of the `true'
% modes. The sign of the modes is also ambiguous.)
pause % Press any key to continue.
echo off
pertrue = 2*pi./abs(angle(EvTrue)); % `true' periods
clc
echo on
% Associated with the eigenmodes are the following oscillation periods:
disp(per)
% The second row contains margins of error for the periods such that
% (per(1,k) +/- per(2,k)) is an approximate 95% confidence interval
% for the period of eigenmode S(:,k). [Note that for a purely
% relaxatory eigenmode, the period is infinite (Inf).] Compare the
% above estimated periods with the `true' periods:
disp(pertrue)
pause % Press any key to get the damping time scales.
echo off
tautrue = -1./log(abs(EvTrue)); % `true' damping time scales
clc
echo on
% The damping times associated with each eigenmode are:
disp(tau)
% with margins of error again in the second row. For comparison, the
% `true' damping times:
disp(tautrue)
pause % Press any key to get the excitation of each eigenmode.
echo off
% Compute `true' excitation of eigenmodes from the designed parameters:
p = 2; % true model order
invStr = inv(Strue); % inverse of matrix with eigenvectors as columns
% covariance matrix of corresponding decoupled AR(1) system
CovDcpld = invStr*[C zeros(2,(p-1)*2); zeros((p-1)*2, p*2)]*invStr';
% diagonal of that covariance matrix
DgCovDcpld = real(diag(CovDcpld))';
% excitation
TrueExctn = DgCovDcpld(1:2*p)./(1-abs(EvTrue).^2);
% normalize excitation
TrueExctn = TrueExctn./sum(TrueExctn);
clc
echo on
% ARMODE has also returned the excitations, measures of the relative
% dynamical importance of the eigenmodes:
disp(exctn)
% Compare the estimated excitations with the `true' excitations
% computed from the parameters used in the simulation:
disp(TrueExctn)
echo off
disp('End')
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
mkdir 'Src'
fi
cd 'Src'
if test -f 'acf.m'
then
echo shar: will not over-write existing file "'acf.m'"
else
cat << "SHAR_EOF" > 'acf.m'
function []=acf(x, k, caption)
%ACF Plot of sample autocorrelation function.
%
% ACF(x) plots the sample autocorrelation function of the univariate
% time series in vector x. By default, the sample autocorrelation
% function is plotted up to lag 25. ACF(x,k) plots the sample
% autocorrelation function up to lag k. ACF(X,k,'name') sets the
% title of the plot to 'name'.
%
% The approximate 95% confidence limits of the autocorrelation
% function of an IID process of the same length as X are also
% displayed. Sample autocorrlations lying outside the 95% confidence
% intervals of an IID process are marked by an asterisk.
%
% ACF requires XCORR from the Signal Processing Toolbox.
%
% See also XCORR.
% Modified 30-Dec-99
% Author: Tapio Schneider
% tapio@cims.nyu.edu
if ~exist('xcorr')
error('ACF requires XCORR from the Signal Processing Toolbox.')
end
[m,n] = size(x);
if (min(m,n) > 1) error('Time series must be univariate.'); end
n = max(m,n);
if (nargin < 3) caption='ACF'; end
if (nargin < 2) k=25; end
if (nargin < 1) error('Usage: ACF(vector).'); end
% Compute autocorrelation sequence with XCORR from the Signal
% Processing Toolbox
cor = xcorr(x,'coeff');
rho = cor(n:n+k); % autocorrelation function up to lag k
abscis = [0:1:k]; % abscissa for plot
% Approximate 95% confidence limits for IID process of the same length
bound = zeros(size(rho));
bound = bound + 1.96/sqrt(n);
% Initialize abscissas and ordinates for ACF values within and outside the
% approximate 95% confidence intervals of IID process
inabsc = zeros(1, k+1);
inl = zeros(1, k+1);
outabsc = zeros(1, k+1);
outl = zeros(1, k+1);
% Find lags within and outside approximate 95% confidence
% intervals; start with lag 0
inl(1) = rho(1); % stores ACF at lags within cfd intervals
inabsc(1) = 0; % lags at which ACF is within cfd intervals
nin = 1; % number of points within confidence intervals
nout = 0; % number of points outside confidence intervals
for i=2:k+1
if abs(rho(i)) > bound(1) % point outside confidence intervals
nout = nout+1;
outl(nout) = rho(i);
outabsc(nout) = i-1;
else % point within confidence intervals
nin = nin+1;
inl(nin) = rho(i);
inabsc(nin) = i-1;
end;
end;
% Plot ACF
plot(abscis, rho, abscis, bound, '-.', abscis, -bound, '-.', ...
outabsc(1:nout), outl(1:nout), '*', ...
inabsc(1:nin), inl(1:nin), 'o');
axis([0 k -1 1])
title(caption)
xlabel('Lag')
ylabel('Autocorrelation function')
SHAR_EOF
fi # end of overwriting check
if test -f 'adjph.m'
then
echo shar: will not over-write existing file "'adjph.m'"
else
cat << "SHAR_EOF" > 'adjph.m'
function ox=adjph(x)
%ADJPH Normalization of columns of a complex matrix.
%
% Given a complex matrix X, OX=ADJPH(X) returns the complex matrix OX
% that is obtained from X by multiplying column vectors of X with
% phase factors exp(i*phi) such that the real part and the imaginary
% part of each column vector of OX are orthogonal and the norm of the
% real part is greater than or equal to the norm of the imaginary
% part.
%
% ADJPH is called by ARMODE.
%
% See also ARMODE.
% Modified 16-Dec-99
% Author: Tapio Schneider
% tapio@cims.nyu.edu
for j = 1:size(x,2)
a = real(x(:,j)); % real part of jth column of x
b = imag(x(:,j)); % imag part of jth column of x
phi = .5*atan( 2*sum(a.*b)/(b'*b-a'*a) );
bnorm = norm(sin(phi).*a+cos(phi).*b); % norm of new imaginary part
anorm = norm(cos(phi).*a-sin(phi).*b); % norm of new real part
if bnorm > anorm
if phi < 0
phi = phi-pi/2;
else
phi = phi+pi/2;
end
end
ox(:,j) = x(:,j).*exp(i*phi);
end
SHAR_EOF
fi # end of overwriting check
if test -f 'arconf.m'
then
echo shar: will not over-write existing file "'arconf.m'"
else
cat << "SHAR_EOF" > 'arconf.m'
function [Aerr, werr]=arconf(A, C, w, th)
%ARCONF Confidence intervals for AR coefficients.
%
% For an AR(p) model that has been fitted with ARFIT,
% [Aerr,werr]=ARCONF(A,C,w,th) computes the margins of error Aerr and
% werr such that (A +/- Aerr) and (w +/- werr) are approximate 95%
% confidence intervals for the elements of the coefficient matrix A
% and for the components of the intercept vector w. The input
% arguments of ARCONF are output of AR.
%
% If no intercept vector w has been fitted with ARFIT (i.e., the flag
% 'zero' was an input argument of ARFIT), then [Aerr]=ARCONF(A,C,th)
% computes the margins of error only for the elements of the
% coefficient matrix A.
%
% The confidence intervals are based on Student's t distribution,
% which for small samples yields only approximate confidence
% intervals. Inferences drawn from small samples must therefore be
% interpreted cautiously.
%
% See also ARFIT.
% Modified 30-Dec-99
% Author: Tapio Schneider
% tapio@cims.nyu.edu
ccoeff = .95; % confidence coefficient
m = size(C,1); % dimension of state space
p = size(A,2)/m; % order of model
if (nargin == 3)
% no intercept vector has been fitted
Aaug = A;
th = w;
w = [];
np = m*p; % number of parameter vectors of size m
else
Aaug = [w A];
np = m*p+1; % number of parameter vectors of size m
end
% number of degrees of freedom for residual covariance matrix
dof = th(1,1);
% quantile of t distribution for given confidence coefficient and dof
t = tquant(dof, .5+ccoeff/2);
% Get matrix Uinv that appears in the covariance matrix of the least squares
% estimator
Uinv = th(2:size(th,1), :);
% Compute approximate confidence intervals for elements of Aaug
Aaug_err = zeros(m, np);
for j=1:m
for k=1:np
Aaug_err(j,k) = t * sqrt( Uinv(k ,k)* C(j,j) );
end
end
if (nargin == 3)
% No intercept vector has been fitted
Aerr = Aaug_err;
else
% An intercept vector has been fitted => return margins of error
% for intercept vector and for AR coefficients separately
werr = Aaug_err(:, 1);
Aerr = Aaug_err(:, 2:np);
end
SHAR_EOF
fi # end of overwriting check
if test -f 'arfit.m'
then
echo shar: will not over-write existing file "'arfit.m'"
else
cat << "SHAR_EOF" > 'arfit.m'
function [w, A, C, sbc, fpe, th]=arfit(v, pmin, pmax, selector, no_const)
%ARFIT Stepwise least squares estimation of multivariate AR model.
%
% [w,A,C,SBC,FPE,th]=ARFIT(v,pmin,pmax) produces estimates of the
% parameters of a multivariate AR model of order p,
%
% v(k,:)' = w' + A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + noise(C),
%
% where p lies between pmin and pmax and is chosen as the optimizer
% of Schwarz's Bayesian Criterion. The input matrix v must contain
% the time series data, with columns of v representing variables
% and rows of v representing observations. ARFIT returns least
% squares estimates of the intercept vector w, of the coefficient
% matrices A1,...,Ap (as A=[A1 ... Ap]), and of the noise covariance
% matrix C.
%
% As order selection criteria, ARFIT computes approximations to
% Schwarz's Bayesian Criterion and to the logarithm of Akaike's Final
% Prediction Error. The order selection criteria for models of order
% pmin:pmax are returned as the vectors SBC and FPE.
%
% The matrix th contains information needed for the computation of
% confidence intervals. ARMODE and ARCONF require th as input
% arguments.
%
% If the optional argument SELECTOR is included in the function call,
% as in ARFIT(v,pmin,pmax,SELECTOR), SELECTOR is used as the order
% selection criterion in determining the optimum model order. The
% three letter string SELECTOR must have one of the two values 'sbc'
% or 'fpe'. (By default, Schwarz's criterion SBC is used.) If the
% bounds pmin and pmax coincide, the order of the estimated model
% is p=pmin=pmax.
%
% If the function call contains the optional argument 'zero' as the
% fourth or fifth argument, a model of the form
%
% v(k,:)' = A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + noise(C)
%
% is fitted to the time series data. That is, the intercept vector w
% is taken to be zero, which amounts to assuming that the AR(p)
% process has zero mean.
% Modified 14-Oct-00
% Authors: Tapio Schneider
% tapio@cims.nyu.edu
%
% Arnold Neumaier
% neum@cma.univie.ac.at
% n: number of observations; m: dimension of state vectors
[n,m] = size(v);
if (pmin ~= round(pmin) | pmax ~= round(pmax))
error('Order must be integer.');
end
if (pmax < pmin)
error('PMAX must be greater than or equal to PMIN.')
end
% set defaults and check for optional arguments
if (nargin == 3) % no optional arguments => set default values
mcor = 1; % fit intercept vector
selector = 'sbc'; % use SBC as order selection criterion
elseif (nargin == 4) % one optional argument
if strcmp(selector, 'zero')
mcor = 0; % no intercept vector to be fitted
selector = 'sbc'; % default order selection
else
mcor = 1; % fit intercept vector
end
elseif (nargin == 5) % two optional arguments
if strcmp(no_const, 'zero')
mcor = 0; % no intercept vector to be fitted
else
error(['Bad argument. Usage: ', ...
'[w,A,C,SBC,FPE,th]=AR(v,pmin,pmax,SELECTOR,''zero'')'])
end
end
ne = n-pmax; % number of block equations of size m
npmax = m*pmax+mcor; % maximum number of parameter vectors of length m
if (ne <= npmax)
error('Time series too short.')
end
% compute QR factorization for model of order pmax
[R, scale] = arqr(v, pmax, mcor);
% compute approximate order selection criteria for models
% of order pmin:pmax
[sbc, fpe] = arord(R, m, mcor, ne, pmin, pmax);
% get index iopt of order that minimizes the order selection
% criterion specified by the variable selector
[val, iopt] = min(eval(selector));
% select order of model
popt = pmin + iopt-1; % estimated optimum order
np = m*popt + mcor; % number of parameter vectors of length m
% decompose R for the optimal model order popt according to
%
% | R11 R12 |
% R=| |
% | 0 R22 |
%
R11 = R(1:np, 1:np);
R12 = R(1:np, npmax+1:npmax+m);
R22 = R(np+1:npmax+m, npmax+1:npmax+m);
% get augmented parameter matrix Aaug=[w A] if mcor=1 and Aaug=A if mcor=0
if (np > 0)
if (mcor == 1)
% improve condition of R11 by re-scaling first column
con = max(scale(2:npmax+m)) / scale(1);
R11(:,1) = R11(:,1)*con;
end;
Aaug = (R11\R12)';
% return coefficient matrix A and intercept vector w separately
if (mcor == 1)
% intercept vector w is first column of Aaug, rest of Aaug is
% coefficient matrix A
w = Aaug(:,1)*con; % undo condition-improving scaling
A = Aaug(:,2:np);
else
% return an intercept vector of zeros
w = zeros(m,1);
A = Aaug;
end
else
% no parameters have been estimated
% => return only covariance matrix estimate and order selection
% criteria for ``zeroth order model''
w = zeros(m,1);
A = [];
end
% return covariance matrix
dof = ne-np; % number of block degrees of freedom
C = R22'*R22./dof; % bias-corrected estimate of covariance matrix
% for later computation of confidence intervals return in th:
% (i) the inverse of U=R11'*R11, which appears in the asymptotic
% covariance matrix of the least squares estimator
% (ii) the number of degrees of freedom of the residual covariance matrix
invR11 = inv(R11);
if (mcor == 1)
% undo condition improving scaling
invR11(1, :) = invR11(1, :) * con;
end
Uinv = invR11*invR11';
th = [dof zeros(1,size(Uinv,2)-1); Uinv];
SHAR_EOF
fi # end of overwriting check
if test -f 'armode.m'
then
echo shar: will not over-write existing file "'armode.m'"
else
cat << "SHAR_EOF" > 'armode.m'
function [S, Serr, per, tau, exctn, lambda] = armode(A, C, th)
%ARMODE Eigendecomposition of AR model.
%
% [S,Serr,per,tau,exctn]=ARMODE(A,C,th) computes the
% eigendecomposition of an AR(p) model that has been fitted using
% ARFIT. The input arguments of ARMODE are output of ARFIT.
%
% The columns of the output matrix S contain the estimated eigenmodes
% of the AR model. The output matrix Serr contains margins of error
% for the components of the estimated eigenmodes S, such that
% (S +/- Serr) are approximate 95% confidence intervals for the
% individual components of the eigenmodes.
%
% The two-row matrices per and tau contain in their first rows the
% estimated oscillation period per(1,k) and the estimated damping
% time tau(1,k) of the eigenmode S(:,k). In their second rows, the
% matrices per and tau contain margins of error for the periods and
% damping times, such that
% ( per(1,k) +/- per(2,k) ) and ( tau(1,k) +/- tau(2,k) )
% are approximate 95% confidence intervals for the period and damping
% time of eigenmode S(:,k).
%
% For a purely relaxatory eigenmode, the period is infinite (Inf).
% For an oscillatory eigenmode, the periods are finite.
%
% The excitation of an eigenmode measures its dynamical importance
% and is returned as a fraction exctn that is normalized such that
% the sum of the excitations of all eigenmodes equals one.
%
% See also ARFIT, ARCONF.
% Modified 13-Oct-00
% Author: Tapio Schneider
% tapio@cims.nyu.edu
ccoeff = .95; % confidence coefficient
m = size(C,1); % dimension of state space
p = size(A,2) / m; % order of model
if p <= 0
error('Order must be greater 0.');
end
% Assemble coefficient matrix of equivalent AR(1) model
A1 = [A; eye((p-1)*m) zeros((p-1)*m,m)];
% Eigenvalues and eigenvectors of coefficient matrix of equivalent
% AR(1) model
[BigS,d] = eig(A1); % columns of BigS are eigenvectors
lambda = diag(d); % vector containing eigenvalues
lambda = lambda(:); % force lambda to be column vector
% Warning if the estimated model is unstable
if any(abs(lambda) > 1)
warning(sprintf(['The estimated AR model is unstable.\n',...
'\t Some excitations may be negative.']))
end
% Fix phase of eigenvectors such that the real part and the
% imaginary part of each vector are orthogonal
BigS = adjph(BigS);
% Return only last m components of each eigenvector
S = BigS((p-1)*m+1:p*m, :);
% Compute inverse of BigS for later use
BigS_inv = inv(BigS);
% Recover the matrix Uinv that appears in the asymptotic covariance
% matrix of the least squares estimator (Uinv is output of AR)
if (size(th,2) == m*p+1)
% The intercept vector has been fitted by AR; in computing
% confidence intervals for the eigenmodes, this vector is
% irrelevant. The first row and first column in Uinv,
% corresponding to elements of the intercept vector, are not
% needed.
Uinv = th(3:size(th,1), 2:size(th,2));
elseif (size(th,2) == m*p)
% No intercept vector has been fitted
Uinv = th(2:size(th,1), :);
else
error('Input arguments of ARMODE must be output of ARFIT.')
end
% Number of degrees of freedom
dof = th(1,1);
% Quantile of t distribution for given confidence coefficient and dof
t = tquant(dof, .5+ccoeff/2);
% Asymptotic covariance matrix of estimator of coefficient matrix A
Sigma_A = kron(Uinv, C);
% Noise covariance matrix of system of relaxators and oscillators
CovDcpld = BigS_inv(:, 1:m) * C * BigS_inv(:, 1:m)';
% For each eigenmode j: compute the period per, the damping time
% tau, and the excitation exctn; also get the margins of error for
% per and tau
for j=1:m*p % eigenmode number
a = real(lambda(j)); % real part of eigenvalue j
b = imag(lambda(j)); % imaginary part of eigenvalue j
abs_lambda_sq= abs(lambda(j))^2; % squared absolute value of eigenvalue j
tau(1,j) = -2/log(abs_lambda_sq);% damping time of eigenmode j
% Excitation of eigenmode j
exctn(j) = real(CovDcpld(j,j) / (1-abs_lambda_sq));
% Assemble derivative of eigenvalue with respect to parameters in
% the coefficient matrix A
dot_lam = zeros(m^2*p, 1);
for k=1:m
dot_lam(k:m:k+(m*p-1)*m) = BigS_inv(j,k) .* BigS(1:m*p,j);
end
dot_a = real(dot_lam); % derivative of real part of lambda(j)
dot_b = imag(dot_lam); % derivative of imag part of lambda(j)
% Derivative of the damping time tau w.r.t. parameters in A
phi = tau(1,j)^2 / abs_lambda_sq * (a*dot_a + b*dot_b);
% Margin of error for damping time tau
tau(2,j) = t * sqrt(phi'*Sigma_A*phi);
% Period of eigenmode j and margin of error for period. (The
% if-statements avoid warning messages that may otherwise result
% from a division by zero)
if (b == 0 & a >= 0) % purely real, nonnegative eigenvalue
per(1,j) = Inf;
per(2,j) = 0;
elseif (b == 0 & a < 0) % purely real, negative eigenvalue
per(1,j) = 2;
per(2,j) = 0;
else % complex eigenvalue
per(1,j) = 2*pi/abs(atan2(b,a));
% Derivative of period with respect to parameters in A
phi = per(1,j)^2 / (2*pi*abs_lambda_sq)*(b*dot_a-a*dot_b);
% Margin of error for period
per(2,j) = t * sqrt(phi'*Sigma_A*phi);
end
end
% Give the excitation as `relative importance' that sums to one
exctn = exctn/sum(exctn);
% Compute confidence intervals for eigenmodes
% -------------------------------------------
% Shorthands for matrix products
XX = real(BigS)'*real(BigS);
YY = imag(BigS)'*imag(BigS);
XY = real(BigS)'*imag(BigS);
% Need confidence intervals only for last m rows of BigS
row1 = (p-1)*m+1; % first row for which confidence interval is needed
mp = m*p; % dimension of equivalent AR(1) model
for k=1:mp % loop over columns of S
for l=row1:mp % loop over rows of S
% evaluate gradient of S_{lk}
for ii=1:m
for jj=1:mp
% compute derivative with respect to A(ii,jj)
zsum = 0;
zkkr = 0; % real part of Z_{kk}
zkki = 0; % imaginary part of Z_{kk}
for j=1:mp
if (j ~= k) % sum up elements that appear in Z_{kk}
zjk = BigS_inv(j,ii)*BigS(jj,k)/(lambda(k)-lambda(j));
zjkr = real(zjk);
zjki = imag(zjk);
zkkr = zkkr+zjki*(XY(k,j)-XY(j,k))-zjkr*(XX(k,j)+YY(k,j));
zkki = zkki+zjki*(YY(k,j)-XX(k,j))-zjkr*(XY(k,j)+XY(j,k));
zsum = zsum+BigS(l,j)*zjk;
end
end
% now add Z_{kk}
zkki = zkki / (XX(k,k)-YY(k,k));
grad_S((jj-1)*m+ii) = zsum+BigS(l,k)*(zkkr+i*zkki);
end
end
Serr(l,k) = t * ( sqrt( real(grad_S)*Sigma_A*real(grad_S)') ...
+ i*sqrt( imag(grad_S)*Sigma_A*imag(grad_S)') );
end
end
% Only return last m*p rows of Serr
Serr = Serr(row1:m*p, :);
SHAR_EOF
fi # end of overwriting check
if test -f 'arord.m'
then
echo shar: will not over-write existing file "'arord.m'"
else
cat << "SHAR_EOF" > 'arord.m'
function [sbc, fpe, logdp, np] = arord(R, m, mcor, ne, pmin, pmax)
%ARORD Evaluates criteria for selecting the order of an AR model.
%
% [SBC,FPE]=ARORD(R,m,mcor,ne,pmin,pmax) returns approximate values
% of the order selection criteria SBC and FPE for models of order
% pmin:pmax. The input matrix R is the upper triangular factor in the
% QR factorization of the AR model; m is the dimension of the state
% vectors; the flag mcor indicates whether or not an intercept vector
% is being fitted; and ne is the number of block equations of size m
% used in the estimation. The returned values of the order selection
% criteria are approximate in that in evaluating a selection
% criterion for an AR model of order p < pmax, pmax-p initial values
% of the given time series are ignored.
%
% ARORD is called by ARFIT.
%
% See also ARFIT, ARQR.
% For testing purposes, ARORD also returns the vectors logdp and np,
% containing the logarithms of the determinants of the (scaled)
% covariance matrix estimates and the number of parameter vectors at
% each order pmin:pmax.
% Modified 17-Dec-99
% Author: Tapio Schneider
% tapio@cims.nyu.edu
imax = pmax-pmin+1; % maximum index of output vectors
% initialize output vectors
sbc = zeros(1, imax); % Schwarz's Bayesian Criterion
fpe = zeros(1, imax); % log of Akaike's Final Prediction Error
logdp = zeros(1, imax); % determinant of (scaled) covariance matrix
np = zeros(1, imax); % number of parameter vectors of length m
np(imax)= m*pmax+mcor;
% Get lower right triangle R22 of R:
%
% | R11 R12 |
% R=| |
% | 0 R22 |
%
R22 = R(np(imax)+1 : np(imax)+m, np(imax)+1 : np(imax)+m);
% From R22, get inverse of residual cross-product matrix for model
% of order pmax
invR22 = inv(R22);
Mp = invR22*invR22';
% For order selection, get determinant of residual cross-product matrix
% logdp = log det(residual cross-product matrix)
logdp(imax) = 2.*log(abs(prod(diag(R22))));
% Compute approximate order selection criteria for models of
% order pmin:pmax
i = imax;
for p = pmax:-1:pmin
np(i) = m*p + mcor; % number of parameter vectors of length m
if p < pmax
% Downdate determinant of residual cross-product matrix
% Rp: Part of R to be added to Cholesky factor of covariance matrix
Rp = R(np(i)+1:np(i)+m, np(imax)+1:np(imax)+m);
% Get Mp, the downdated inverse of the residual cross-product
% matrix, using the Woodbury formula
L = chol(eye(m) + Rp*Mp*Rp')';
N = L \ Rp*Mp;
Mp = Mp - N'*N;
% Get downdated logarithm of determinant
logdp(i) = logdp(i+1) + 2.* log(abs(prod(diag(L))));
end
% Schwarz's Bayesian Criterion
sbc(i) = logdp(i)/m - log(ne) * (ne-np(i))/ne;
% logarithm of Akaike's Final Prediction Error
fpe(i) = logdp(i)/m - log(ne*(ne-np(i))/(ne+np(i)));
% Modified Schwarz criterion (MSC):
% msc(i) = logdp(i)/m - (log(ne) - 2.5) * (1 - 2.5*np(i)/(ne-np(i)));
i = i-1; % go to next lower order
end
SHAR_EOF
fi # end of overwriting check
if test -f 'arqr.m'
then
echo shar: will not over-write existing file "'arqr.m'"
else
cat << "SHAR_EOF" > 'arqr.m'
function [R, scale]=arqr(v, p, mcor)
%ARQR QR factorization for least squares estimation of AR model.
%
% [R, SCALE]=ARQR(v,p,mcor) computes the QR factorization needed in
% the least squares estimation of parameters of an AR(p) model. If
% the input flag mcor equals one, a vector of intercept terms is
% being fitted. If mcor equals zero, the process v is assumed to have
% mean zero. The output argument R is the upper triangular matrix
% appearing in the QR factorization of the AR model, and SCALE is a
% vector of scaling factors used to regularize the QR factorization.
%
% ARQR is called by ARFIT.
%
% See also ARFIT.
% Modified 29-Dec-99
% Author: Tapio Schneider
% tapio@cims.nyu.edu
% n: number of time steps; m: dimension of state vectors
[n,m] = size(v);
ne = n-p; % number of block equations of size m
np = m*p+mcor; % number of parameter vectors of size m
% If the intercept vector w is to be fitted, least squares (LS)
% estimation proceeds by solving the normal equations for the linear
% regression model
%
% v(k,:)' = Aaug*u(k,:)' + noise(C)
%
% with Aaug=[w A] and `predictors'
%
% u(k,:) = [1 v(k-1,:) ... v(k-p,:)].
%
% If the process mean is taken to be zero, the augmented coefficient
% matrix is Aaug=A, and the regression model
%
% u(k,:) = [v(k-1,:) ... v(k-p,:)]
%
% is fitted.
% The number np is the dimension of the `predictors' u(k).
% Assemble the data matrix K (of which a QR factorization will be computed)
K = zeros(ne,np+m); % initialize K
if (mcor == 1)
% first column of K consists of ones for estimation of intercept vector w
K(:,1) = ones(ne,1);
end
% Assemble `predictors' u in K
for j=1:p
K(:, mcor+m*(j-1)+1:mcor+m*j) = [v(p-j+1:n-j, :)];
end
% Add `observations' v (left hand side of regression model) to K
K(:,np+1:np+m) = [v(p+1:n, :)];
% Compute regularized QR factorization of K: The regularization
% parameter delta is chosen according to Higham's (1996) Theorem
% 10.7 on the stability of a Cholesky factorization. Replace the
% regularization parameter delta below by a parameter that depends
% on the observational error if the observational error dominates
% the rounding error (cf. Neumaier, A. and Schneider, T.,
% "Estimation of parameters and eigenmodes of multivariate
% autoregressive models", ACM Trans. Math. Softw., 2001).
q = np + m; % number of columns of K
delta = (q^2 + q + 1)*eps; % Higham's choice for a Cholesky factorization
scale = sqrt(delta)*sqrt(sum(K.^2));
R = triu(qr([K; diag(scale)]));
SHAR_EOF
fi # end of overwriting check
if test -f 'arres.m'
then
echo shar: will not over-write existing file "'arres.m'"
else
cat << "SHAR_EOF" > 'arres.m'
function [siglev,res]=arres(w,A,v,k)
%ARRES Test of residuals of fitted AR model.
%
% [siglev,res]=ARRES(w,A,v) computes the time series of residuals
%
% res(k,:)' = v(k+p,:)'- w - A1*v(k+p-1,:)' - ... - Ap*v(k,:)'
%
% of an AR(p) model with A=[A1 ... Ap].
%
% Also returned is the significance level siglev of the modified
% Li-McLeod portmanteau (LMP) statistic.
%
% Correlation matrices for the LMP statistic are computed up to lag
% k=20, which can be changed to lag k by using
% [siglev,res]=ARRES(w,A,v,k).
% Modified 17-Dec-99
% Author: Tapio Schneider
% tapio@cims.nyu.edu
%
% Reference:
% Li, W. K., and A. I. McLeod, 1981: Distribution of the
% Residual Autocorrelations in Multivariate ARMA Time Series
% Models, J. Roy. Stat. Soc. B, 43, 231--239.
m = size(v,2); % dimension of state vectors
p = size(A,2)/m; % order of model
n = length(v); % number of observations
nres = n-p; % number of residuals
% Default value for k
if (nargin < 4)
k = 20;
end
if (k <= p) % check if k is in valid range
error('Maximum lag of residual correlation matrices too small.');
end
if (k >= nres)
error('Maximum lag of residual correlation matrices too large.');
end
w = w(:)'; % force w to be row vector
% Get time series of residuals
l = 1:nres; % vectorized loop l=1,...,nres
res(l,:) = v(l+p,:) - ones(nres,1)*w;
for j=1:p
res(l,:) = res(l,:) - v(l-j+p,:)*A(:, (j-1)*m+1:j*m)';
end
% end of loop over l
% Center residuals by subtraction of the mean
res = res - ones(nres,1)*mean(res);
% Compute lag zero correlation matrix of the residuals
c0 = res'*res;
d = diag(c0);
dd = sqrt(d*d');
c0 = c0./dd;
% Get "covariance matrix" in LMP statistic
c0_inv= inv(c0); % inverse of lag 0 correlation matrix
rr = kron(c0_inv, c0_inv); % "covariance matrix" in LMP statistic
% Initialize LMP statistic and correlation matrices
lmp = 0; % LMP statistic
cl = zeros(m,m); % correlation matrix
x = zeros(m*m,1); % correlation matrix arranged as vector
% Compute modified Li-McLeod portmanteau statistic
for l=1:k
cl = (res(1:nres-l, :)'*res(l+1:nres,:))./dd; % lag l correlation matrix
x = reshape(cl,m*m, 1); % arrange cl as vector by stacking columns
lmp = lmp + x'*rr*x; % sum up LMP statistic
end
lmp = n*lmp + m^2*k*(k+1)/2/n; % add remaining term and scale
dof_lmp = m^2*(k-p); % degrees of freedom for LMP statistic
% Significance level with which hypothesis of uncorrelatedness is rejected
siglev = 1 - gammainc(lmp/2, dof_lmp/2);
SHAR_EOF
fi # end of overwriting check
if test -f 'arsim.m'
then
echo shar: will not over-write existing file "'arsim.m'"
else
cat << "SHAR_EOF" > 'arsim.m'
function [v]=arsim(w,A,C,n,ndisc)
%ARSIM Simulation of AR process.
%
% v=ARSIM(w,A,C,n) simulates n time steps of the AR(p) process
%
% v(k,:)' = w' + A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + eta(k,:)',
%
% where A=[A1 ... Ap] is the coefficient matrix, and w is a vector of
% intercept terms that is included to allow for a nonzero mean of the
% process. The vectors eta(k,:) are independent Gaussian noise
% vectors with mean zero and covariance matrix C.
%
% The p vectors of initial values for the simulation are taken to
% be equal to the mean value of the process. (The process mean is
% calculated from the parameters A and w.) To avoid spin-up effects,
% the first 10^3 time steps are discarded. Alternatively,
% ARSIM(w,A,C,n,ndisc) discards the first ndisc time steps.
% Modified 13-Oct-00
% Author: Tapio Schneider
% tapio@cims.nyu.edu
m = size(C,1); % dimension of state vectors
p = size(A,2)/m; % order of process
if (p ~= round(p))
error('Bad arguments.');
end
if (length(w) ~= m | min(size(w)) ~= 1)
error('Dimensions of arguments are mutually incompatible.')
end
w = w(:)'; % force w to be row vector
% Check whether specified model is stable
A1 = [A; eye((p-1)*m) zeros((p-1)*m,m)];
lambda = eig(A1);
if any(abs(lambda) > 1)
warning('The specified AR model is unstable.')
end
% Discard the first ndisc time steps; if ndisc is not given as input
% argument, use default
if (nargin < 5)
ndisc = 10^3;
end
% Compute Cholesky factor of covariance matrix C
[R, err]= chol(C); % R is upper triangular
if err ~= 0
error('Covariance matrix not positive definite.')
end
% Get ndisc+n independent Gaussian pseudo-random vectors with
% covariance matrix C=R'*R
randvec = randn([ndisc+n,m])*R;
% Add intercept vector to random vectors
randvec = randvec + ones(ndisc+n,1)*w;
% Get transpose of system matrix A (use transpose in simulation because
% we want to obtain the states as row vectors)
AT = A';
% Take the p initial values of the simulation to equal the process mean,
% which is calculated from the parameters A and w
if any(w)
% Process has nonzero mean mval = inv(B)*w' where
% B = eye(m) - A1 -... - Ap;
% Assemble B
B = eye(m);
for j=1:p
B = B - A(:, (j-1)*m+1:j*m);
end
% Get mean value of process
mval = w / B';
% The optimal forecast of the next state given the p previous
% states is stored in the vector x. The vector x is initialized
% with the process mean.
x = ones(p,1)*mval;
else
% Process has zero mean
x = zeros(p,m);
end
% Initialize state vectors
u = [x; zeros(ndisc+n,m)];
% Simulate n+ndisc observations. In order to be able to make use of
% Matlab's vectorization capabilities, the cases p=1 and p>1 must be
% treated separately.
if p==1
for k=2:ndisc+n+1;
x(1,:) = u(k-1,:)*AT;
u(k,:) = x + randvec(k-1,:);
end
else
for k=p+1:ndisc+n+p;
for j=1:p;
x(j,:) = u(k-j,:)*AT((j-1)*m+1:j*m,:);
end
u(k,:) = sum(x)+randvec(k-p,:);
end
end
% return only the last n simulated state vectors
v = u(ndisc+p+1:ndisc+n+p,:);
SHAR_EOF
fi # end of overwriting check
if test -f 'tquant.m'
then
echo shar: will not over-write existing file "'tquant.m'"
else
cat << "SHAR_EOF" > 'tquant.m'
function t=tquant(n, p)
%TQUANT Quantiles of Student's t distribution
%
% TQUANT(n, p) is the p-quantile of a t distributed random variable
% with n degrees of freedom; that is, TQUANT(n, p) is the value below
% which 100p percent of the t distribution with n degrees of freedom
% lies.
% Modified 15-Apr-00
% Author: Tapio Schneider
% tapio@cims.nyu.edu
% References:
% L. Devroye, 1986: "Non-Uniform Random Variate Generation", Springer
%
% M. Abramowitz and I. A. Stegun, 1964: "Handbook of Mathematical
% Functions"
%
% See also: tcdf.m in the Matlab Statistics Toolbox (evaluates
% cumulative distribution function of Student's t)
if (n ~= round(n) | n < 1)
error('Usage: TQUANT(n,p) - Degrees of freedom n must be positive integer.')
end
if (p<0 | p>1)
error('Usage: TQUANT(n,p) - Probability p must be in [0,1].')
elseif p == 1
t = Inf;
return
elseif p == 0
t = -Inf;
return
end
if n == 1
% Cauchy distribution (cf. Devroye [1986, pp. 29 and 450])
t = tan(pi*(p-.5));
elseif p >= 0.5
% positive t-values (cf. M. Abramowitz and I. A. Stegun [1964,
% Chapter 26])
b0 = [0, 1];
f = inline('1 - betainc(b, n/2, .5)/2 - p', 'b', 'n', 'p');
%opt = optimset('Display', 'off'); % does not work on Windows/Mac
%b = fzero(f, b0, opt, n, p);
% old calling sequence (for Windows/Mac compatibility):
b = fzero(f, b0, eps, 0, n, p);
t = sqrt(n/b-n);
else
% negative t-values
b0 = [0, 1];
f = inline('betainc(b, n/2, .5)/2 - p', 'b', 'n', 'p');
%opt = optimset('Display', 'off'); % does not work on Windows/Mac
%b = fzero(f, b0, opt, n, p);
% old calling sequence (for Windows/Mac compatibility):
b = fzero(f, b0, eps, 0, n, p);
t = -sqrt(n/b-n);
end
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
# End of shell archive
exit 0