hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_tsa_uni_arima_estim (g13ae)

Purpose

nag_tsa_uni_arima_estim (g13ae) fits a seasonal autoregressive integrated moving average (ARIMA) model to an observed time series, using a nonlinear least squares procedure incorporating backforecasting. Parameter estimates are obtained, together with appropriate standard errors. The residual series is returned, and information for use in forecasting the time series is produced for use by the functions nag_tsa_uni_arima_update (g13ag) and nag_tsa_uni_arima_forecast_state (g13ah).
The estimation procedure is iterative, starting with initial parameter values such as may be obtained using nag_tsa_uni_arima_prelim (g13ad). It continues until a specified convergence criterion is satisfied, or until a specified number of iterations has been carried out. The progress of the procedure can be monitored by means of a user-supplied function.

Syntax

[par, c, icount, ex, exr, al, s, g, sd, h, st, nst, itc, zsp, isf, ifail] = g13ae(mr, par, c, x, iex, igh, ist, piv, kpiv, zsp, kzsp, 'npar', npar, 'kfc', kfc, 'nx', nx, 'nit', nit)
[par, c, icount, ex, exr, al, s, g, sd, h, st, nst, itc, zsp, isf, ifail] = nag_tsa_uni_arima_estim(mr, par, c, x, iex, igh, ist, piv, kpiv, zsp, kzsp, 'npar', npar, 'kfc', kfc, 'nx', nx, 'nit', nit)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 23: nit and kfc optional
.

Description

The time series x1,x2,,xnx1,x2,,xn supplied to nag_tsa_uni_arima_estim (g13ae) is assumed to follow a seasonal autoregressive integrated moving average (ARIMA) model defined as follows:
dsDxtc = wt,
dsDxt-c=wt,
where dsDxtdsDxt is the result of applying non-seasonal differencing of order dd and seasonal differencing of seasonality ss and order DD to the series xtxt, as outlined in the description of nag_tsa_uni_diff (g13aa). The differenced series is then of length N = ndN=n-d, where d = d + (D × s)d=d+(D×s) is the generalized order of differencing. The scalar cc is the expected value of the differenced series, and the series w1,w2,,wNw1,w2,,wN follows a zero-mean stationary autoregressive moving average (ARMA) model defined by a pair of recurrence equations. These express wtwt in terms of an uncorrelated series atat, via an intermediate series etet. The first equation describes the seasonal structure:
wt = Φ1wts + Φ2wt2 × s + + ΦPwtP × s + etΘ1etsΘ2et2 × sΘQetQ × s.
wt=Φ1wt-s+Φ2wt-2×s++ΦPwt-P×s+et-Θ1et-s-Θ2et-2×s--ΘQet-Q×s.
The second equation describes the non-seasonal structure. If the model is purely non-seasonal the first equation is redundant and etet above is equated with wtwt:
et = φ1et1 + φ2et2 + + φpetp + atθ1at1θ2at2θqatq.
et=ϕ1et-1+ϕ2et-2++ϕpet-p+at-θ1at-1-θ2at-2--θqat-q.
Estimates of the model parameters defined by
φ1,φ2,,φp,θ1,θ2,,θq,
Φ1,Φ2,,ΦP,Θ1,Θ2,,ΘQ
ϕ1,ϕ2,,ϕp,θ1,θ2,,θq, Φ1,Φ2,,ΦP,Θ1,Θ2,,ΘQ
and (optionally) cc are obtained by minimizing a quadratic form in the vector w = (w1,w2,,wN)w=(w1,w2,,wN).
This is QF = wV1wQF=wV-1w, where VV is the covariance matrix of ww, and is a function of the model parameters. This matrix is not explicitly evaluated, since QFQF may be expressed as a ‘sum of squares’ function. When moving average parameters θiθi or ΘiΘi are present, so that the generalized moving average order q = q + s × Qq=q+s×Q is positive, backforecasts w1q,w2q,,w0w1-q,w2-q,,w0 are introduced as nuisance parameters. The ‘sum of squares’ function may then be written as
N q
S(pm) = at2bt2,
t = 1q t = 1qp
S(pm)=t=1-qNat2-t=1-q-p -qbt2,
where pmpm is a combined vector of parameters, consisting of the backforecasts followed by the ARMA model parameters.
The terms atat correspond to the ARMA model residual series atat, and p = p + s × Pp=p+s×P is the generalized autoregressive order. The terms btbt are only present if autoregressive parameters are in the model, and serve to correct for transient errors introduced at the start of the autoregression.
The equations defining atat and btbt are precisely:
For all four of these equations, the following conditions hold:
Minimization of SS with respect to pmpm uses an extension of the algorithm of Marquardt (1963).
The first derivatives of SS with respect to the parameters are calculated as
2 × at × at,i2bt × bt,i = 2 × Gi,
2×at×at,i-2bt×bt,i=2×Gi,
where at,iat,i and bt,ibt,i are derivatives of atat and btbt with respect to the iith parameter.
The second derivative of SS is approximated by
2 × at,i × at,j2 × bt,i × bt,j = 2 × Hij.
2×at,i×at,j-2×bt,i×bt,j=2×Hij.
Successive parameter iterates are obtained by calculating a vector of corrections dpmdpm by solving the equations
(H + α × D) × dpm = G,
(H+α×D)×dpm=-G,
where GG is a vector with elements GiGi, HH is a matrix with elements HijHij, αα is a scalar used to control the search and DD is the diagonal matrix of HH.
The new parameter values are then pm + dpmpm+dpm.
The scalar αα controls the step size, to which it is inversely related.
If a step results in new parameter values which give a reduced value of SS, then αα is reduced by a factor ββ. If a step results in new parameter values which give an increased value of SS, or in ARMA model parameters which in any way contravene the stationarity and invertibility conditions, then the new parameters are rejected, αα is increased by the factor ββ, and the revised equations are solved for a new parameter correction.
This action is repeated until either a reduced value of SS is obtained, or αα reaches the limit of 109109, which is used to indicate a failure of the search procedure.
This failure may be due to a badly conditioned sum of squares function or to too strict a convergence criterion. Convergence is deemed to have occurred if the fractional reduction in the residual sum of squares in successive iterations is less than a value γγ, while α < 1.0α<1.0.
The stationarity and invertibility conditions are tested to within a specified tolerance multiple δδ of machine accuracy. Upon convergence, or completion of the specified maximum number of iterations without convergence, statistical properties of the estimates are derived. In the latter case the sequence of iterates should be checked to ensure that convergence is adequate for practical purposes, otherwise these properties are not reliable.
The estimated residual variance is
erv = Smin / df ,
erv = Smin / df ,
where SminSmin is the final value of SS, and the residual number of degrees of freedom is given by
df = NpqPQ (1​ if ​c​ is estimated) .
df = N-p-q-P-Q ( -1 ​ if ​ c ​ is estimated ) .
The covariance matrix of the vector of estimates pmpm is given by
erv × H1 ,
erv × H-1 ,
where HH is evaluated at the final parameter values.
From this expression are derived the vector of standard deviations, and the correlation matrix for the whole parameter set. These are asymptotic approximations.
The differenced series wtwt (now uncorrected for the constant), intermediate series etet and residual series atat are all available upon completion of the iterations over the range (extended by backforecasts)
t = 1q,2q,,N.
t=1-q,2-q,,N.
The values atat can only properly be interpreted as residuals for t1 + pqt1+p-q, as the earlier values are corrupted by transients if p > 0p>0.
In consequence of the manner in which differencing is implemented, the residual atat is the one step ahead forecast error for xt + dxt+d.
For convenient application in forecasting, the following quantities constitute the ‘state set’, which contains the minimum amount of time series information needed to construct forecasts:
(i) the differenced series wtwt, for (Ns × P) < tN(N-s×P)<tN,
(ii) the dd values required to reconstitute the original series xtxt from the differenced series wtwt,
(iii) the intermediate series etet, for (Nmax (p, Q × s )) < t N ( N - max(p, Q × s ) ) < t N ,
(iv) the residual series atat, for (Nq) < tN(N-q)<tN.
This state set is available upon completion of the iterations. The function may be used purely for the construction of this state set, given a previously estimated model and time series xtxt, by requesting zero iterations. Backforecasts are estimated, but the model parameter values are unchanged. If later observations become available and it is desired to update the state set, nag_tsa_uni_arima_update (g13ag) can be used.

References

Box G E P and Jenkins G M (1976) Time Series Analysis: Forecasting and Control (Revised Edition) Holden–Day
Marquardt D W (1963) An algorithm for least squares estimation of nonlinear parameters J. Soc. Indust. Appl. Math. 11 431

Parameters

Compulsory Input Parameters

1:     mr(77) – int64int32nag_int array
The orders vector (p,d,q,P,D,Q,s)(p,d,q,P,D,Q,s) of the ARIMA model whose parameters are to be estimated. pp, qq, PP and QQ refer respectively to the number of autoregressive (φϕ), moving average (θ)(θ), seasonal autoregressive (ΦΦ) and seasonal moving average (ΘΘ) parameters. dd, DD and ss refer respectively to the order of non-seasonal differencing, the order of seasonal differencing and the seasonal period.
Constraints:
  • pp, dd, qq, PP, DD, QQ, s0s0;
  • p + q + P + Q > 0p+q+P+Q>0;
  • s1s1;
  • if s = 0s=0, P + D + Q = 0P+D+Q=0;
  • if s > 1s>1, P + D + Q > 0P+D+Q>0;
  • d + s × (P + D)nd+s×(P+D)n;
  • p + dq + s × (P + DQ)np+d-q+s×(P+D-Q)n.
2:     par(npar) – double array
npar, the dimension of the array, must satisfy the constraint npar = p + q + P + Qnpar=p+q+P+Q.
The initial estimates of the pp values of the φϕ parameters, the qq values of the θθ parameters, the PP values of the ΦΦ parameters and the QQ values of the ΘΘ parameters, in that order.
3:     c – double scalar
If kfc = 0kfc=0, c must contain the expected value, cc, of the differenced series.
If kfc = 1kfc=1, c must contain an initial estimate of cc.
4:     x(nx) – double array
The nn values of the original undifferenced time series.
5:     iex – int64int32nag_int scalar
The dimension of the arrays ex, exr and al as declared in the (sub)program from which nag_tsa_uni_arima_estim (g13ae) is called.
Constraint: iexq + (Q × s) + niexq+(Q×s)+n, which is equivalent to the exit value of icount(4)icount4.
6:     igh – int64int32nag_int scalar
The dimension of the arrays g and sd and the second dimension of the arrays h and hc as declared in the (sub)program from which nag_tsa_uni_arima_estim (g13ae) is called.
Constraint: ighq + (Q × s) + npar + kfcighq+(Q×s)+npar+kfc which is equivalent to the exit value of icount(6)icount6.
7:     ist – int64int32nag_int scalar
The dimension of the array st as declared in the (sub)program from which nag_tsa_uni_arima_estim (g13ae) is called.
Constraint: ist(P × s) + d + (D × s) + q + max (p,Q × s)ist(P×s)+d+(D×s)+q+max(p,Q×s).
8:     piv – function handle or string containing name of m-file
piv is used to monitor the progress of the optimization.
piv(mr, par, npar, c, kfc, icount, s, g, h, ldh, igh, itc, zsp)

Input Parameters

1:     mr(77) – int64int32nag_int array
2:     par(npar) – double array
3:     npar – int64int32nag_int scalar
4:     c – double scalar
5:     kfc – int64int32nag_int scalar
6:     icount(66) – int64int32nag_int array
7:     s – double scalar
8:     g(igh) – double array
9:     h(ldh,igh) – double array
10:   ldh – int64int32nag_int scalar
11:   igh – int64int32nag_int scalar
12:   itc – int64int32nag_int scalar
13:   zsp(44) – double array
All the parameters are defined as for nag_tsa_uni_arima_estim (g13ae) itself.
If kpiv = 0kpiv=0 a dummy piv must be supplied.
9:     kpiv – int64int32nag_int scalar
Must be nonzero if the progress of the optimization is to be monitored using piv. Otherwise kpiv must contain 00.
10:   zsp(44) – double array
When kzsp = 1kzsp=1, the first four elements of zsp must contain the four values used to guide the search procedure. These are as follows.
zsp(1)zsp1 contains αα, the value used to constrain the magnitude of the search procedure steps.
zsp(2)zsp2 contains ββ, the multiplier which regulates the value αα.
zsp(3)zsp3 contains δδ, the value of the stationarity and invertibility test tolerance factor.
zsp(4)zsp4 contains γγ, the value of the convergence criterion.
If kzsp1kzsp1 on entry, default values for zsp are supplied by the function.
These are 0.0010.001, 10.010.0, 1000.01000.0 and max (100 × machine precision,0.0000001)max(100×machine precision, 0.0000001) respectively.
Constraint: if kzsp = 1kzsp=1, zsp(1) > 0.0zsp1>0.0, zsp(2) > 1.0zsp2>1.0, zsp(3)1.0zsp31.0, 0zsp(4) < 1.00zsp4<1.0.
11:   kzsp – int64int32nag_int scalar
The value 11 if the function is to use the input values of zsp, and any other value if the default values of zsp are to be used.

Optional Input Parameters

1:     npar – int64int32nag_int scalar
Default: The dimension of the array par.
The total number of φϕ, θθ, ΦΦ and ΘΘ parameters to be estimated.
Constraint: npar = p + q + P + Qnpar=p+q+P+Q.
2:     kfc – int64int32nag_int scalar
Must be set to 11 if the constant, cc, is to be estimated and 00 if it is to be held fixed at its initial value.
Default: 11 
Constraint: kfc = 0kfc=0 or 11.
3:     nx – int64int32nag_int scalar
Default: The dimension of the array x.
nn, the length of the original undifferenced time series.
4:     nit – int64int32nag_int scalar
The maximum number of iterations to be performed.
Default: 100100 
Constraint: nit0nit0.

Input Parameters Omitted from the MATLAB Interface

ldh wa iwa hc

Output Parameters

1:     par(npar) – double array
The latest values of the estimates of these parameters.
2:     c – double scalar
If kfc = 0kfc=0, c is unchanged.
If kfc = 1kfc=1, c contains the latest estimate of cc.
Therefore, if c and kfc are both zero on entry, there is no constant correction.
3:     icount(66) – int64int32nag_int array
Size of various output arrays.
icount(1)icount1
Contains q + (Q × s)q+(Q×s), the number of backforecasts.
icount(2)icount2
Contains nd(D × s)n-d-(D×s), the number of differenced values.
icount(3)icount3
Contains d + (D × s)d+(D×s), the number of values of reconstitution information.
icount(4)icount4
Contains n + q + (Q × s)n+q+(Q×s), the number of values held in each of the series ex, exr and al.
icount(5)icount5
Contains nd(D × s)pqPQkfcn-d-(D×s)-p-q-P-Q-kfc, the number of degrees of freedom associated with SS.
icount(6)icount6
Contains icount(1) + npar + kfcicount1+npar+kfc, the number of parameters being estimated.
These values are always computed regardless of the exit value of ifail.
4:     ex(iex) – double array
The extended differenced series which is made up of:
icount(1)icount1 backforecast values of the differenced series.
icount(2)icount2 actual values of the differenced series.
icount(3)icount3 values of reconstitution information.
The total number of these values held in ex is icount(4)icount4.
If the function exits because of a faulty input parameter, the contents of ex will be indeterminate.
5:     exr(iex) – double array
The values of the model residuals which is made up of:
icount(1)icount1 residuals corresponding to the backforecasts in the differenced series.
icount(2)icount2 residuals corresponding to the actual values in the differenced series.
The remaining icount(3)icount3 values contain zeros.
If the function exits with ifail holding a value other than 00 or 99, the contents of exr will be indeterminate.
6:     al(iex) – double array
The intermediate series which is made up of:
icount(1)icount1 intermediate series values corresponding to the backforecasts in the differenced series.
icount(2)icount2 intermediate series values corresponding to the actual values in the differenced series.
The remaining icount(3)icount3 values contain zeros.
If the function exits with ifail0ifail0, the contents of al will be indeterminate.
7:     s – double scalar
The residual sum of squares after the latest series of parameter estimates has been incorporated into the model. If the function exits with a faulty input parameter, s contains zero.
8:     g(igh) – double array
The latest value of the derivatives of SS with respect to each of the parameters being estimated (backforecasts, par parameters, and where relevant the constant – in that order). The contents of g will be indeterminate if the function exits with a faulty input parameter.
9:     sd(igh) – double array
The standard deviations corresponding to each of the parameters being estimated (backforecasts, par parameters, and where relevant the constant, in that order).
If the function exits with ifail containing a value other than 00 or 99, or if the required number of iterations is zero, the contents of sd will be indeterminate.
10:   h(ldh,igh) – double array
ldh1 + q + (q × s) + npar + kfcldh1+q+(q×s)+npar+kfc, which is equivalent to the exit value of icount(6)icount6.
The second derivative of SS and correlation coefficients.
(a) the latest values of an approximation to the second derivative of SS with respect to each of the (q + Q × s + npar + kfc)(q+Q×s+npar+kfc) parameters being estimated (backforecasts, par parameters, and where relevant the constant – in that order), and
(b) the correlation coefficients relating to each pair of these parameters.
These are held in a matrix defined by the first (q + Q × s + npar + kfc)(q+Q×s+npar+kfc) rows and the first (q + Q × s + npar + kfc)(q+Q×s+npar+kfc) columns of h. (Note that icount(6)icount6 contains the value of this expression.) The values of (a) are contained in the upper triangle, and the values of (b) in the strictly lower triangle.
These correlation coefficients are zero during intermediate printout using piv, and indeterminate if ifail contains on exit a value other than 00 or 99.
All the contents of h are indeterminate if the required number of iterations are zero. The (q + (Q × s) + npar + kfc + 1)(q+(Q×s)+npar+kfc+1)th row of h is used internally as workspace.
11:   st(ist) – double array
The nst values of the state set array. If the function exits with ifail containing a value other than 00 or 99, the contents of st will be indeterminate.
12:   nst – int64int32nag_int scalar
The number of values in the state set array st.
13:   itc – int64int32nag_int scalar
The number of iterations performed.
14:   zsp(44) – double array
zsp contains the values, default or otherwise, used by the function.
15:   isf(44) – int64int32nag_int array
Contains success/failure indicators, one for each of the four types of parameter in the model (autoregressive, moving average, seasonal autoregressive, seasonal moving average), in that order.
Each indicator has the interpretation:
2-2 On entry parameters of this type have initial estimates which do not satisfy the stationarity or invertibility test conditions.
1-1 The search procedure has failed to converge because the latest set of parameter estimates of this type is invalid.
0-0 No parameter of this type is in the model.
1-1 Valid final estimates for parameters of this type have been obtained.
16:   ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_tsa_uni_arima_estim (g13ae) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

  ifail = 1ifail=1
On entry,nparp + q + P + Qnparp+q+P+Q,
orthe orders vector mr is invalid (check it against the constraints in Section [Parameters]),
orkfc0kfc0 or 11.
  ifail = 2ifail=2
On entry, nxdD × snpar + kfcnx-d-D×snpar+kfc, i.e., the number of terms in the differenced series is not greater than the number of parameters in the model. The model is over-parameterised.
  ifail = 3ifail=3
On entry,one or more of the user-supplied criteria for controlling the iterative process are invalid,
ornit < 0nit<0,
orif kzsp = 1kzsp=1, zsp(1)0.0zsp10.0;
orif kzsp = 1kzsp=1, zsp(2)1.0zsp21.0;
orif kzsp = 1kzsp=1, zsp(3) < 1.0zsp3<1.0;
orif kzsp = 1kzsp=1, zsp(4) < 0.0zsp4<0.0;
orif kzsp = 1kzsp=1, zsp(4)1.0zsp41.0.
  ifail = 4ifail=4
On entry, the state set array st is too small. The output value of nst contains the required value (see the description of ist in Section [Parameters] for the formula).
  ifail = 5ifail=5
On entry, the workspace array wa is too small. Check the value of iwa against the constraints in Section [Parameters].
W ifail = 6ifail=6
On entry,iex < q + (Q × s) + nxiex<q+(Q×s)+nx,
origh < q + (Q × s) + npar + kfcigh<q+(Q×s)+npar+kfc,
orldhq + (Q × s) + npar + kfcldhq+(Q×s)+npar+kfc.
W ifail = 7ifail=7
This indicates a failure in the search procedure, with zsp(1)1.0e09zsp11.0e09.
Some output parameters may contain meaningful values; see Section [Parameters] for details.
W ifail = 8ifail=8
This indicates a failure to invert HH.
Some output parameters may contain meaningful values; see Section [Parameters] for details.
W ifail = 9ifail=9
This indicates a failure in nag_linsys_real_posdef_solve_1rhs (f04as) which is used to solve the equations giving the latest estimates of the backforecasts.
Some output parameters may contain meaningful values; see Section [Parameters] for details.
W ifail = 10ifail=10
Satisfactory parameter estimates could not be obtained for all parameter types in the model. Inspect array isf for further information on the parameter type(s) in error.

Accuracy

The computations are believed to be stable.

Further Comments

The time taken by nag_tsa_uni_arima_estim (g13ae) is approximately proportional to nx × itc × (q + Q × s + npar + kfc)2nx×itc× (q+Q×s+npar+kfc) 2.

Example

function nag_tsa_uni_arima_estim_example
mr = [int64(1);1;2;0;0;0;0];
par = [0; 0; 0];
c = 0;
x = [-217;
     -177;
     -166;
     -136;
     -110;
     -95;
     -64;
     -37;
     -14;
     -25;
     -51;
     -62;
     -73;
     -88;
     -113;
     -120;
     -83;
     -33;
     -19;
     21;
     17;
     44;
     44;
     78;
     88;
     122;
     126;
     114;
     85;
     64];
iex = int64(32);
igh = int64(6);
ist = int64(4);
kpiv = int64(0);
zsp = [0.001; 10; 1000; 0.0001];
kzsp = int64(1);
[parOut, cOut, icount, ex, exr, al, s, g, sd, h, st, nst, itc, zspOut, isf, ifail] = ...
    nag_tsa_uni_arima_estim(mr, par, c, x, iex, igh, ist, @piv, kpiv, zsp, kzsp)

function [] = piv(mr, par, npar, c, kfc, icount, s, g, h, ih, igh, itc, zsp)

  fprintf('Iteration %d  residual sum f squares = %16.4', itc, s);
 

parOut =

   -0.0547
   -0.5568
   -0.6636


cOut =

    9.9807


icount =

                    2
                   29
                    1
                   32
                   25
                    6


ex =

   19.5250
    5.8753
   40.0000
   11.0000
   30.0000
   26.0000
   15.0000
   31.0000
   27.0000
   23.0000
  -11.0000
  -26.0000
  -11.0000
  -11.0000
  -15.0000
  -25.0000
   -7.0000
   37.0000
   50.0000
   14.0000
   40.0000
   -4.0000
   27.0000
         0
   34.0000
   10.0000
   34.0000
    4.0000
  -12.0000
  -29.0000
  -21.0000
   64.0000


exr =

   19.5250
   -3.9279
   19.5711
   -5.6291
   10.2221
   15.1582
   -9.3276
   16.4285
   15.2115
   -5.4211
  -27.3444
  -18.3061
    5.3890
  -12.9812
  -22.4767
  -15.2183
    4.4944
   33.6867
   19.7586
  -27.1470
   32.2426
  -12.2765
    1.6941
   -1.8465
   23.3772
  -10.4576
   14.3302
   -5.7061
  -28.6401
  -20.4502
   -2.7215
         0


al =

   19.5250
    5.8753
   30.0193
    1.0193
   20.0193
   16.0193
    5.0193
   21.0193
   17.0193
   13.0193
  -20.9807
  -35.9807
  -20.9807
  -20.9807
  -24.9807
  -34.9807
  -16.9807
   27.0193
   40.0193
    4.0193
   30.0193
  -13.9807
   17.0193
   -9.9807
   24.0193
    0.0193
   24.0193
   -5.9807
  -21.9807
  -38.9807
  -30.9807
         0


s =

   9.3979e+03


g =

   -0.1512
   -0.2343
   -6.4097
   13.5617
  -72.6232
   -0.1642


sd =

   14.8379
   15.1887
    0.3507
    0.2709
    0.1695
    7.3893


h =

   1.0e+04 *

    0.0002   -0.0001    0.0000    0.0002   -0.0001    0.0000
    0.0000    0.0002   -0.0000   -0.0000    0.0002    0.0001
   -0.0000    0.0000    0.9042   -0.9682    0.0546    0.0001
   -0.0000    0.0000    0.0001    1.7031   -0.5676    0.0007
   -0.0000   -0.0000    0.0000    0.0000    1.7028    0.0006
   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0007
   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0000


st =

   64.0000
  -30.9807
  -20.4502
   -2.7215


nst =

                    4


itc =

                   16


zspOut =

   1.0e+03 *

    0.0000
    0.0100
    1.0000
    0.0000


isf =

                    1
                    1
                    0
                    0


ifail =

                    0


function g13ae_example
mr = [int64(1);1;2;0;0;0;0];
par = [0; 0; 0];
c = 0;
x = [-217;
     -177;
     -166;
     -136;
     -110;
     -95;
     -64;
     -37;
     -14;
     -25;
     -51;
     -62;
     -73;
     -88;
     -113;
     -120;
     -83;
     -33;
     -19;
     21;
     17;
     44;
     44;
     78;
     88;
     122;
     126;
     114;
     85;
     64];
iex = int64(32);
igh = int64(6);
ist = int64(4);
kpiv = int64(0);
zsp = [0.001; 10; 1000; 0.0001];
kzsp = int64(1);
[parOut, cOut, icount, ex, exr, al, s, g, sd, h, st, nst, itc, zspOut, isf, ifail] = ...
    g13ae(mr, par, c, x, iex, igh, ist, @piv, kpiv, zsp, kzsp)

function [] = piv(mr, par, npar, c, kfc, icount, s, g, h, ih, igh, itc, zsp)

  fprintf('Iteration %d  residual sum f squares = %16.4', itc, s);
 

parOut =

   -0.0547
   -0.5568
   -0.6636


cOut =

    9.9807


icount =

                    2
                   29
                    1
                   32
                   25
                    6


ex =

   19.5250
    5.8753
   40.0000
   11.0000
   30.0000
   26.0000
   15.0000
   31.0000
   27.0000
   23.0000
  -11.0000
  -26.0000
  -11.0000
  -11.0000
  -15.0000
  -25.0000
   -7.0000
   37.0000
   50.0000
   14.0000
   40.0000
   -4.0000
   27.0000
         0
   34.0000
   10.0000
   34.0000
    4.0000
  -12.0000
  -29.0000
  -21.0000
   64.0000


exr =

   19.5250
   -3.9279
   19.5711
   -5.6291
   10.2221
   15.1582
   -9.3276
   16.4285
   15.2115
   -5.4211
  -27.3444
  -18.3061
    5.3890
  -12.9812
  -22.4767
  -15.2183
    4.4944
   33.6867
   19.7586
  -27.1470
   32.2426
  -12.2765
    1.6941
   -1.8465
   23.3772
  -10.4576
   14.3302
   -5.7061
  -28.6401
  -20.4502
   -2.7215
         0


al =

   19.5250
    5.8753
   30.0193
    1.0193
   20.0193
   16.0193
    5.0193
   21.0193
   17.0193
   13.0193
  -20.9807
  -35.9807
  -20.9807
  -20.9807
  -24.9807
  -34.9807
  -16.9807
   27.0193
   40.0193
    4.0193
   30.0193
  -13.9807
   17.0193
   -9.9807
   24.0193
    0.0193
   24.0193
   -5.9807
  -21.9807
  -38.9807
  -30.9807
         0


s =

   9.3979e+03


g =

   -0.1512
   -0.2343
   -6.4097
   13.5617
  -72.6232
   -0.1642


sd =

   14.8379
   15.1887
    0.3507
    0.2709
    0.1695
    7.3893


h =

   1.0e+04 *

    0.0002   -0.0001    0.0000    0.0002   -0.0001    0.0000
    0.0000    0.0002   -0.0000   -0.0000    0.0002    0.0001
   -0.0000    0.0000    0.9042   -0.9682    0.0546    0.0001
   -0.0000    0.0000    0.0001    1.7031   -0.5676    0.0007
   -0.0000   -0.0000    0.0000    0.0000    1.7028    0.0006
   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0007
   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0000


st =

   64.0000
  -30.9807
  -20.4502
   -2.7215


nst =

                    4


itc =

                   16


zspOut =

   1.0e+03 *

    0.0000
    0.0100
    1.0000
    0.0000


isf =

                    1
                    1
                    0
                    0


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013