AUTO_PARM
Estimates structural breaks in non‑stationary univariate time series.
Required Arguments
Y — Array containing the time series. (Input)
NPCS — Number of requested/estimated pieces or segments of the time series. NPCS is considered input when IFITONLY = 1. (Input/Output)
ARP — A pointer to an array of size NPCS × 2. ARP is considered input when IFITONLY = 1. (Input/Output)
Column Index | Description |
|---|---|
1 | Requested/estimated break points |
2 | AR order for each segment |
ARFIT — A pointer to an array of size NPCS × MAXARORDER containing the AR coefficient estimates for each segment. ARFIT (i,j) is the j‑th coefficient for segment i where i=1, NPCS and j=l, ARP(i,2). (Output)
Note that the intercept is not reported.
Note that the intercept is not reported.
ARSTAT — A pointer to an array of size NPCS × 2. (Output)
Column Index | Description |
|---|---|
1 | Likelihood values for each of the fitted AR models |
2 | Residual variances for each of the fitted AR models |
SCVAL — Final value of the selection criterion. (Output)
Optional
MAXARORDER—Maximum order to consider for each AR model. (Input)
Default: MAXARORDER = 20.
Default: MAXARORDER = 20.
IMTH — Method of estimation. (Input)
Value | Method |
|---|---|
0 | Yule‑Walker |
1 | Least Squares |
2 | Burg |
Default: IMTH = 0.
ISELCRI — Selection criterion. (Input)
Value | Method |
|---|---|
0 | Minimum Description Length (MDL) |
1 | Akaike’s Information Criterion (AIC) |
Default: ISELCRI = 0.
ILIKE — Likelihood computation method. (Input)
Value | Method |
|---|---|
0 | Exact |
1 | Approximate |
Default: ILIKE = 0.
IFITONLY — Option to only fit the specified model. (Input)
Value | Action |
|---|---|
0 | No, do all the computations |
1 | Yes, only fit the specified model |
Default: IFITONLY = 0.
IPRINT — Printing option. (Input)
Value | Action |
|---|---|
0 | No printing |
1 | Prints final results only |
2 | Prints intermediate and final results |
Default: IPRINT = 0.
ISEED — Seed of the random number generator. (Input)
For the same data and parameter settings, AUTO_PARM will return the same results each time if ISEED = 1, 2 … 2147483646. If ISEED = 0, the system clock will be used to generate a seed. The result will be nondeterministic.
Default: ISEED = 0.
For the same data and parameter settings, AUTO_PARM will return the same results each time if ISEED = 1, 2 … 2147483646. If ISEED = 0, the system clock will be used to generate a seed. The result will be nondeterministic.
Default: ISEED = 0.
Note: The following arguments are for setting up and running the embedded Genetic Algorithm. In most situations, the default values should be used for these arguments. Users may wish to change some or all for testing or research purposes. |
PDISTN — Array of length MAXARORDER + 1 giving the probability distribution over the AR order variable p = 0,…, MAXARORDER. (Input)
j = 1,…, MAXARORDER + 1 is used to randomly assign an AR order to breakpoint position j for a given chromosome. PDISTN(j) > = 0 and if SUM(PDISTN) is not equal to 1, the values will be normalized, i.e., PDISTN(j) = PDISTN(j)/SUM(PDISTN).
Default: PDISTN (j) = 1/(MAXARORDER + 1) for all j.
j = 1,…, MAXARORDER + 1 is used to randomly assign an AR order to breakpoint position j for a given chromosome. PDISTN(j) > = 0 and if SUM(PDISTN) is not equal to 1, the values will be normalized, i.e., PDISTN(j) = PDISTN(j)/SUM(PDISTN).
Default: PDISTN (j) = 1/(MAXARORDER + 1) for all j.
MSPAN — Array of length MAXARORDER + 1 containing minimum number of observations required for valid estimates of AR model with order p = 0, …, MAXARORDER. (Input)
Default: MSPAN (p+ 1) = 2 * (number of parameters) + 2 = 2 * (p + 2) + 2.
Default: MSPAN (p+ 1) = 2 * (number of parameters) + 2 = 2 * (p + 2) + 2.
GAPARM — Array of length 4 containing parameters that control the behavior of the genetic algorithm. These values should be strictly greater than zero and less than one to avoid unexpected results. (Input)
GAPARM (1) — Probability used to set initial break points in a chromosome.
Default: MIN (MSPAN) / size(Y).
Default: MIN (MSPAN) / size(Y).
GAPARM (2) — Probability of crossover used to decide between a crossover and a mutation.
Default: 1 – MIN (MSPAN) / size(Y).
Default: 1 – MIN (MSPAN) / size(Y).
GAPARM (3) — In the mutation operation, probability an AR(p) model is enforced at the current position.
Default: 0.4.
Default: 0.4.
GAPARM (4) — In the mutation operation, probability a break point is disallowed at the current position.
Default: 0.3.
Default: 0.3.
Note: GAPARM(3) and GAPARM(4) must be valid probabilities and their sum must be between 0 and 1. 1 – GAPARM(3)– GAPARM(4) is the probability that the chromosome j inherits the parent’s chromosome j. |
ISLAND — Array of length 5 containing the migration policy parameters. (Input)
ISLAND (1) — Number of islands.
Default: 40.
Default: 40.
ISLAND (2) — Number of generations that pass before migration occurs. Note that the convergence of the algorithm is revised whether migrations take place or not (see argument ISLAND(5)).
Default: 5.
Default: 5.
ISLAND (3) — Number of subjects that migrate at each migration event.
Default: 2.
Default: 2.
ISLAND (4) — Population size (number of chromosomes) per island.
Default: 40.
Default: 40.
ISLAND (5) = Migration flag. If 1, migration is performed.
Default: 1.
Default: 1.
MAXMIG — Maximum number of times that migrations may take place before the routine is stopped if convergence has not occurred. (Input)
Default: MAXMIG = 20.
Default: MAXMIG = 20.
STOPITERS — Number of iterations. The routine will declare convergence and stop the iterations if the criterion value (MDL/AIC) has not changed after STOPITERS consecutive migrations. Otherwise, the algorithm will declare non‑convergence and stop after MAXMIG migrations have taken place. See also MAXMIG and ISLAND(2). Note that logically, STOPITERS < MAXMIG. (Input)
Default: STOPITERS = 10.
Default: STOPITERS = 10.
Interface
Generic: CALL AUTO_PARM (Y, NPCS, ARP, ARFIT, ARSTAT, SCVAL [, …])
Specific: The specific interface names are S_AUTO_PARM and D_AUTO_PARM.
Description
Routine AUTO_PARM estimates the structural breaks of a non‑stationary time series using, with permission from the Authors, the method developed by Davis et al (2006). AUTO_PARM estimates a partition of the time index and models the time series in each segment as a separate auto‑regressive (AR(p)) process. The routine returns the estimated breakpoints, the estimated AR(p) models, and supporting statistics.
For the observed time series , the problem is to find m(m+1 = NPCS), the number of breaks, their locations, (ARP(:,1)), and Pj(ARP(:,2)), j = 1, …, m + 1, the order of the AR process in which the j-th segment is modeled. That is, for (for convenience, and ) where {Xt,j} is an AR process of order Pj.
and ɛt, the noise sequence, is i.i.d. with mean 0 and variance 1. Note that a series with m breaks will have m + 1 segments.
The vector completely specifies a piecewise AR model. To estimate this vector AUTO_PARM optimizes, with respect to this vector, one of two selection critieria: the first is a Minimum Description Length (MDL) criterion, and the second is the Akaike’s Information Criterion (AIC). The MDL is defined as
while the AIC criterion is given by
where, given a candidate value of the vector , L is the likelihood of the fitted piecewise AR model evaluated at the parameter estimates,
The parameters of the j-th AR segment are estimated by the choice of one of three estimation methods: Yule-Walker, Burg, or Least Squares.
For simplicity, assume the mean of each series is 0 and that the errors are Gaussian. Then, the piecewise AR model has Gaussian likelihood
where is the variance-covariance of the j-th AR segment (of order ) and is the vector of observations of the j-th segment, i.e., 
To find the minimizer {m, j, pj}* of either MDL or AIC, AUTO_PARM employs a Genetic Algorithm with islands, migration, cross‑over and mutations. See Davis et.al. (2006) for further details.
Comments
1. AUTO_PARM approximates locally stationary time series by independent auto‑regressive processes. Experimental results suggest that AUTO_PARM gives reasonable estimates of the structural breaks of a given time series, even if the segment series are not autoregressive. Also, based on experimental results, MDL gives better results than AIC as a selection criterion.
2. Informational Error
Type | Code | Description |
|---|---|---|
3 | 1 | ISEED has been set out of range. ISEED is being reset to 123457. |
3 | 2 | MAXMIG migrations were reached in the genetic algorithm before the selection criterion value converged. Try increasing MAXMIG or using the double precision routine. |
Examples
The examples below illustrate different scenarios using AUTO_PARM. The example series used in each case is the airline demand data (Box, Jenkins and Reinsel, 1994), which gives monthly total demand for the period January 1949 through December 1960. Each scenario sets the optional argument, ISEED = 123457.
use gdata_int
use auto_parm_int
use diff_int
use umach_int
implicit none
! Specifications for local variables
! arguments for auto_parm
integer :: npcs, maxarorder, ifitonly, iprint, island(5), iseed
integer, pointer :: arp(:,:)
real(kind(1e0)) :: x(144,1), scval
real(kind(1e0)), pointer, dimension(:,:) :: arfit, arstat
! Other locals
integer :: n, nvar, idata, nlost, iper(1), iord(1), nout
real(kind(1e0)) :: dx(144,1)
call umach(2, nout)
! Get the data
iseed = 123457
idata = 4
nullify(arp)
nullify(arfit)
nullify(arstat)
call gdata(idata, x, n, nvar)
! Example 1: Use defaults and print
! final results
write(nout,*)'Example 1: Use defaults'
write(nout,*)
call auto_parm(x(:,1), npcs, arp, arfit, arstat, scval, &
iseed=iseed, iprint=1)
! Example 2: Differenced series
! set period for the difference.
! iper is in years for this data set
write(nout,*)
write(nout,*)'Example 2: Differenced series'
write(nout,*)
! Set the order for the difference.
iper(1) = 1
iord(1) = 1
! Get differenced series dx
call diff(x(:,1),iper,iord,n,dx(:,1),nlost=nlost)
! Compare results on the differenced
! series
deallocate(arp)
deallocate(arfit)
deallocate(arstat)
call auto_parm(dx((1+nlost):n,1), npcs, arp, arfit, arstat, &
scval, iseed=iseed, iprint=1)
! Example 3: Original series
! lower maximum AR order
write(nout,*)
write(nout,*)'Example 3: Original series, lower order allowed'
write(nout,*)
maxarorder = 5
deallocate(arp)
deallocate(arfit)
deallocate(arstat)
call auto_parm(x(:,1), npcs, arp, arfit, arstat, scval, &
maxarorder=maxarorder, iseed=iseed, iprint=1)
! Example 4: differenced series, lower
! maximum AR order
write(nout,*)
write(nout,*)'Example 4: Differenced series, lower maximum'// &
' AR order'
write(nout,*)
deallocate(arp)
deallocate(arfit)
deallocate(arstat)
call auto_parm(dx((1+nlost):n,1), npcs, arp, arfit, arstat, &
scval, maxarorder=maxarorder, iseed=iseed, &
iprint=1)
! Example 5: Original series, force
! fit the segments
! Fit the specified model only at the
! break points
write(nout,*)
write(nout,*)'Example 5: Original series'
write(nout,*)'Force fit the segments at the break points'
write(nout,*)
npcs = 2
deallocate(arp)
allocate(arp(npcs,2))
arp(1,2) = 2
arp(2,2) = 1
arp(1,1) = 1
arp(2,1) = 60
deallocate(arfit)
deallocate(arstat)
call auto_parm(x(:,1), npcs, arp, arfit, arstat, scval, &
ifitonly=1, iseed=iseed, iprint=1)
end
Output
Example 1: Use defaults
============== final results ===============
number of pieces: 2
selection criteria value: 684.242
total time: 0.7198125 conv: 1
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 1 0.77542
186.945 355.025
ARP( 2,1) ARP( 2,2) ARFIT( 2,:) ARSTAT( 2,1) ARSTAT( 2,2)
44 13 1.03701
-0.07802
-0.03890
-0.03452
0.11961
-0.12852
0.01990
-0.04886
0.08090
-0.13118
0.22122
0.53863
-0.61515
486.665 691.471
Example 2: Differenced series
============== final results ===============
number of pieces: 1
selection criteria value: 624.284
total time: 0.7204375 conv: 1
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 12 -0.02842
-0.22436
-0.16846
-0.24267
-0.10573
-0.22429
-0.12126
-0.26446
-0.07087
-0.24327
-0.07136
0.57129
619.321 297.351
Example 3: Original series, lower order allowed
============== final results ===============
number of pieces: 2
selection criteria value: 705.297
total time: 0.3136875 conv: 1
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 1 0.89533
270.393 333.564
ARP( 2,1) ARP( 2,2) ARFIT( 2,:) ARSTAT( 2,1) ARSTAT( 2,2)
63 2 1.19788
-0.35922
424.270 1632.337
Example 4: Differenced series, lower maximum AR order
============== final results ===============
number of pieces: 2
selection criteria value: 698.359
total time: 0.2981875 conv: 1
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 0 -------
335.565 357.388
ARP( 2,1) ARP( 2,2) ARFIT( 2,:) ARSTAT( 2,1) ARSTAT( 2,2)
77 1 0.33310
352.175 1786.345
Example 5: Original series
Force fit the segments at the break points
============== final results ===============
number of pieces: 2
selection criteria value: 712.521
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 2 1.12156
-0.24876
258.192 313.889
ARP( 2,1) ARP( 2,2) ARFIT( 2,:) ARSTAT( 2,1) ARSTAT( 2,2)
60 1 0.88605
443.696 1937.635
