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_spectrum_lag (g13ca)

▸▿ Contents

    1  Purpose
    2  Syntax
    3  Description
    4  References
▸▿ 5  Parameters
    6  Error Indicators and Warnings
    7  Accuracy
    8  Further Comments
    9  Example

Purpose

nag_tsa_uni_spectrum_lag (g13ca) calculates the smoothed sample spectrum of a univariate time series using one of four lag windows – rectangular, Bartlett, Tukey or Parzen window.

Syntax

[c, xg, ng, stats, ifail] = g13ca(nx, mtx, px, iw, mw, ic, c, kc, l, lg, xg, 'nc', nc, 'nxg', nxg)
[c, xg, ng, stats, ifail] = nag_tsa_uni_spectrum_lag(nx, mtx, px, iw, mw, ic, c, kc, l, lg, xg, 'nc', nc, 'nxg', nxg)

Description

The smoothed sample spectrum is defined as
f^ω=12π C0+2∑k=1 M-1wkCkcosωk ,  
where M is the window width, and is calculated for frequency values
ωi=2πiL,  i=0,1,…,L/2,  
where  denotes the integer part.
The autocovariances Ck may be supplied by you, or constructed from a time series x1,x2,…,xn, as
Ck=1n∑t=1 n-kxtxt+k,  
the fast Fourier transform (FFT) being used to carry out the convolution in this formula.
The time series may be mean or trend corrected (by classical least squares), and tapered before calculation of the covariances, the tapering factors being those of the split cosine bell:
121-cosπ t-12/T, 1≤t≤T 121-cosπ n-t+12/T, n+ 1-T≤t≤n 1, otherwise,  
where T= np2  and p is the tapering proportion.
The smoothing window is defined by
wk=W kM ,  k≤M-1,  
which for the various windows is defined over 0≤α<1 by
rectangular:
Wα=1  
Bartlett:
Wα = 1-α  
Tukey:
Wα=121+cosπα  
Parzen:
Wα= 1- 6α2+ 6α3, 0≤α≤12 Wα= 2 1-α 3, 12<α< 1.  
The sampling distribution of f^ω is approximately that of a scaled χd2 variate, whose degrees of freedom d is provided by the function, together with multiplying limits mu, ml from which approximate 95% confidence intervals for the true spectrum fω may be constructed as ml × f ^ ω , mu × f ^ ω . Alternatively, log f^ω may be returned, with additive limits.
The bandwidth b of the corresponding smoothing window in the frequency domain is also provided. Spectrum estimates separated by (angular) frequencies much greater than b may be assumed to be independent.

References

Bloomfield P (1976) Fourier Analysis of Time Series: An Introduction Wiley
Jenkins G M and Watts D G (1968) Spectral Analysis and its Applications Holden–Day

Parameters

Compulsory Input Parameters

1:     nx – int64int32nag_int scalar
n, the length of the time series.
Constraint: nx≥1.
2:     mtx – int64int32nag_int scalar
If covariances are to be calculated by the function (ic=0), mtx must specify whether the data are to be initially mean or trend corrected.
mtx=0
For no correction.
mtx=1
For mean correction.
mtx=2
For trend correction.
Constraint: if ic=0, 0≤mtx≤2
If covariances are supplied (ic≠0), mtx is not used.
3:     px – double scalar
If covariances are to be calculated by the function (ic=0), px must specify the proportion of the data (totalled over both ends) to be initially tapered by the split cosine bell taper.
If covariances are supplied ic≠0, px must specify the proportion of data tapered before the supplied covariances were calculated and after any mean or trend correction. px is required for the calculation of output statistics. A value of 0.0 implies no tapering.
Constraint: 0.0≤px≤1.0.
4:     iw – int64int32nag_int scalar
The choice of lag window.
iw=1
Rectangular.
iw=2
Bartlett.
iw=3
Tukey.
iw=4
Parzen.
Constraint: 1≤iw≤4.
5:     mw – int64int32nag_int scalar
M, the ‘cut-off’ point of the lag window. Windowed covariances at lag M or greater are zero.
Constraint: 1≤mw≤nx.
6:     ic – int64int32nag_int scalar
Indicates whether covariances are to be calculated in the function or supplied in the call to the function.
ic=0
Covariances are to be calculated.
ic≠0
Covariances are to be supplied.
7:     cnc – double array
If ic≠0, c must contain the nc covariances for lags from 0 to nc-1, otherwise c need not be set.
8:     kc – int64int32nag_int scalar
If ic=0, kc must specify the order of the fast Fourier transform (FFT) used to calculate the covariances. kc should be a product of small primes such as 2m where m is the smallest integer such that 2m≥nx+nc, provided m≤20.
If ic≠0, that is covariances are supplied, kc is not used.
Constraint: kc≥nx+nc. The largest prime factor of kc must not exceed 19, and the total number of prime factors of kc, counting repetitions, must not exceed 20. These two restrictions are imposed by the internal FFT algorithm used.
9:     l – int64int32nag_int scalar
L, the frequency division of the spectral estimates as 2πL . Therefore it is also the order of the FFT used to construct the sample spectrum from the covariances. l should be a product of small primes such as 2m where m is the smallest integer such that 2m≥2M-1, provided m≤20.
Constraint: l≥2×mw-1. The largest prime factor of l must not exceed 19, and the total number of prime factors of l, counting repetitions, must not exceed 20. These two restrictions are imposed by the internal FFT algorithm used.
10:   lg – int64int32nag_int scalar
Indicates whether unlogged or logged spectral estimates and confidence limits are required.
lg=0
Unlogged.
lg≠0
Logged.
11:   xgnxg – double array
If the covariances are to be calculated, then xg must contain the nx data points. If covariances are supplied, xg may contain any values.

Optional Input Parameters

1:     nc – int64int32nag_int scalar
Default: the dimension of the array c.
The number of covariances to be calculated in the function or supplied in the call to the function.
Constraint: mw≤nc≤nx.
2:     nxg – int64int32nag_int scalar
Default: the dimension of the array xg.
The dimension of the array xg.
Constraints:
  • if ic=0, nxg≥maxkc,l;
  • if ic≠0, nxg≥l.

Output Parameters

1:     cnc – double array
If ic=0, c will contain the nc calculated covariances.
If ic≠0, the contents of c will be unchanged.
2:     xgnxg – double array
Contains the ng spectral estimates, f^ωi, for i=0,1,…,L/2 in xg1 to xgng respectively (logged if lg=1). The elements xgi, for i=ng+1,…,nxg contain 0.0.
3:     ng – int64int32nag_int scalar
The number of spectral estimates, L/2+1, in xg.
4:     stats4 – double array
Four associated statistics. These are the degrees of freedom in stats1, the lower and upper 95% confidence limit factors in stats2 and stats3 respectively (logged if lg=1), and the bandwidth in stats4.
5:     ifail – int64int32nag_int scalar
ifail=0 unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and 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=1
On entry,nx<1,
ormtx<0 and ic=0,
ormtx>2 and ic=0,
orpx<0.0,
orpx>1.0,
oriw<1,
oriw>4,
ormw<1,
ormw>nx,
ornc<mw,
ornc>nx,
ornxg<maxkc,l and ic=0,
ornxg<l and ic≠0.
   ifail=2
On entry,kc<nx+nc,
orkc has a prime factor exceeding 19,
orkc has more than 20 prime factors, counting repetitions.
This error only occurs when ic=0.
   ifail=3
On entry,l<2×mw-1,
orl has a prime factor exceeding 19,
orl has more than 20 prime factors, counting repetitions.
W  ifail=4
One or more spectral estimates are negative. Unlogged spectral estimates are returned in xg, and the degrees of freedom, unlogged confidence limit factors and bandwidth in stats.
W  ifail=5
The calculation of confidence limit factors has failed. This error will not normally occur. Spectral estimates (logged if requested) are returned in xg, and degrees of freedom and bandwidth in stats.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   ifail=-999
Dynamic memory allocation failed.

Accuracy

The FFT is a numerically stable process, and any errors introduced during the computation will normally be insignificant compared with uncertainty in the data.

Further Comments

nag_tsa_uni_spectrum_lag (g13ca) carries out two FFTs of length kc to calculate the covariances and one FFT of length l to calculate the sample spectrum. The time taken by the function for an FFT of length n is approximately proportional to nlogn (but see Further Comments in nag_sum_fft_realherm_1d (c06pa) for further details).

Example

This example reads a time series of length 256. It selects the mean correction option, a tapering proportion of 0.1, the Parzen smoothing window and a cut-off point for the window at lag 100. It chooses to have 100 auto-covariances calculated and unlogged spectral estimates at a frequency division of 2π/200. It then calls nag_tsa_uni_spectrum_lag (g13ca) to calculate the univariate spectrum and statistics and prints the autocovariances and the spectrum together with its 95% confidence multiplying limits.

Open in the MATLAB editor: g13ca_example

function g13ca_example


fprintf('g13ca example results\n\n');

% Sizes
nx  = int64(256);
nc  = 100;
kc  = int64(360);

%Smoothing parameters
mtx = int64(1);
ic  = int64(0);
px  = 0.1;
iw  = int64(4);
mw  = int64(100);
l   = int64(200);
lg  = int64(0);


% Arrays
c  = zeros(nc, 1);
xg = zeros(kc,1);
xg(1:nx) = ...
     [ ...
        5.0  11.0  16.0  23.0  36.0  58.0  29.0  20.0  10.0   8.0 ...
        3.0   0.0   0.0   2.0  11.0  27.0  47.0  63.0  60.0  39.0 ... 
       28.0  26.0  22.0  11.0  21.0  40.0  78.0 122.0 103.0  73.0 ...
       47.0  35.0  11.0   5.0  16.0  34.0  70.0  81.0 111.0 101.0 ...
       73.0  40.0  20.0  16.0   5.0  11.0  22.0  40.0  60.0  80.9 ...
       83.4  47.7  47.8  30.7  12.2   9.6  10.2  32.4  47.6  54.0 ...
       62.9  85.9  61.2  45.1  36.4  20.9  11.4  37.8  69.8 106.1 ...
      100.8  81.6  66.5  34.8  30.6   7.0  19.8  92.5 154.4 125.9 ...
       84.8  68.1  38.5  22.8  10.2  24.1  82.9 132.0 130.9 118.1 ...
       89.9  66.6  60.0  46.9  41.0  21.3  16.0   6.4   4.1   6.8 ...
       14.5  34.0  45.0  43.1  47.5  42.2  28.1  10.1   8.1   2.5 ...
        0.0   1.4   5.0  12.2  13.9  35.4  45.8  41.1  30.1  23.9 ...
       15.6   6.6   4.0   1.8   8.5  16.6  36.3  49.6  64.2  67.0 ...
       70.9  47.8  27.5   8.5  13.2  56.9 121.5 138.3 103.2  85.7 ...
       64.6  36.7  24.2  10.7  15.0  40.1  61.5  98.5 124.7  96.3 ...
       66.6  64.5  54.1  39.0  20.6   6.7   4.3  22.7  54.8  93.8 ...
       95.8  77.2  59.1  44.0  47.0  30.5  16.3   7.3  37.6  74.0 ...
      139.0 111.2 101.6  66.2  44.7  17.0  11.3  12.4   3.4   6.0 ...
       32.3  54.3  59.7  63.7  63.5  52.2  25.4  13.1   6.8   6.3 ...
        7.1  35.6  73.0  85.1  78.0  64.0  41.8  26.2  26.7  12.1 ...
        9.5   2.7   5.0  24.4  42.0  63.5  53.8  62.0  48.5  43.9 ...
       18.6   5.7   3.6   1.4   9.6  47.4  57.1 103.9  80.6  63.6 ...
       37.6  26.1  14.2   5.8  16.7  44.3  63.9  69.0  77.8  64.9 ...
       35.7  21.2  11.1   5.7   8.7  36.1  79.7 114.4 109.6  88.8 ...
       67.8  47.5  30.6  16.3   9.6  33.2  92.6 151.6 136.3 134.7 ...
       83.9  69.4  31.5  13.9   4.4  38.0];

% Calculate smoothed spectrum
[c, xg, ng, stats, ifail] = ...
g13ca( ...
       nx, mtx, px, iw, mw, ic, c, kc, l, lg, xg);

% Display results
disp('Covariances');
for j = 1:6:nc
  fprintf('%11.4f', c(j:min(j+5,nc)));
  fprintf('\n');
end
fprintf('\nDegrees of freedom = %4.1f      Bandwidth = %7.4f\n\n', ...
        stats(1), stats(4));
fprintf('95 percent confidence limits -  Lower = %7.4f  Upper = %7.4f\n', ...
        stats(2:3));
fprintf('\n      Spectrum       Spectrum       Spectrum       Spectrum\n');
fprintf('      estimate       estimate       estimate       estimate\n');
result = [double([1:ng]); xg(1:ng)'];
for j = 1:4:ng
  fprintf('%4d%10.4f', result(:,j:min(j+3,ng)));
  fprintf('\n');
end



g13ca example results

Covariances
  1152.9733   937.3289   494.9243    14.8648  -342.8548  -514.6479
  -469.2733  -236.6896   109.0608   441.3498   637.4571   641.9954
   454.0505   154.5960  -136.8016  -343.3911  -421.8441  -374.4095
  -241.1943   -55.6140   129.4067   267.4248   311.8293   230.2807
    56.4402  -146.4689  -320.9948  -406.4077  -375.6384  -273.5936
  -132.6214    11.0791   126.4843   171.3391   122.6284   -11.5482
  -169.2623  -285.2358  -331.4567  -302.2945  -215.4832  -107.8732
    -3.4126    73.2521    98.0831    71.8949    17.0985   -27.5632
   -76.7900  -110.5354  -126.1383  -121.1043  -103.9362   -67.4619
   -10.8678    58.5009   116.4587   140.0961   129.5928    66.3211
   -35.5487  -135.3894  -203.7149  -216.2161  -152.7723   -30.4361
    99.3397   188.9594   204.9047   148.4056    34.4975  -103.7840
  -208.5982  -252.4128  -223.7600  -120.8640    23.3565   156.0956
   227.7642   228.5123   172.3820    87.4911   -21.2170  -117.5282
  -176.3634  -165.1218   -75.1308    67.1634   195.7290   279.3039
   290.8258   225.3811   104.0784   -44.4731  -162.7355  -207.7480
  -165.2444   -48.5473   118.8872   265.0045

Degrees of freedom =  9.0      Bandwidth =  0.1165

95 percent confidence limits -  Lower =  0.4731  Upper =  3.3329

      Spectrum       Spectrum       Spectrum       Spectrum
      estimate       estimate       estimate       estimate
   1  210.4696   2  428.2020   3  810.1419   4  922.5900
   5  706.1605   6  393.4052   7  207.6481   8  179.0657
   9  170.1320  10  133.0442  11  103.6752  12  103.0644
  13  141.5173  14  194.3041  15  266.5730  16  437.0181
  17  985.3130  18 2023.1574  19 2681.8980  20 2363.7439
  21 1669.9001  22 1012.1320  23  561.4822  24  467.2741
  25  441.9977  26  300.1985  27  172.0184  28  114.7823
  29   79.1533  30   49.4882  31   27.0902  32   16.8081
  33   27.5111  34   59.4429  35   97.0145  36  119.3664
  37  116.6737  38   87.3142  39   54.9570  40   42.9781
  41   46.6097  42   53.6206  43   50.6050  44   36.7780
  45   25.6285  46   24.8555  47   30.2626  48   31.5642
  49   27.3351  50   22.4443  51   18.5418  52   15.2425
  53   12.0207  54   12.6846  55   18.3975  56   19.3058
  57   12.6103  58    7.9511  59    7.1333  60    5.4996
  61    3.4182  62    3.2359  63    5.3836  64    8.5225
  65   10.0610  66    7.9483  67    4.2261  68    3.2631
  69    5.5751  70    7.8491  71    9.3694  72   11.0791
  73   10.1386  74    6.3158  75    3.6375  76    2.6561
  77    1.8026  78    1.0103  79    1.0693  80    2.3950
  81    4.0822  82    4.6221  83    4.0672  84    3.8460
  85    4.8489  86    6.3964  87    6.4762  88    4.9457
  89    4.4444  90    5.2131  91    5.0389  92    4.6141
  93    5.8722  94    7.9268  95    7.9486  96    5.7854
  97    4.5495  98    5.2696  99    6.3893 100    6.5216
 101    6.2129
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–2015