Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_tsa_uni_spectrum_lag (g13ca)

## 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 $\left[\right]$ denotes the integer part.
The autocovariances ${C}_{k}$ may be supplied by you, or constructed from a time series ${x}_{1},{x}_{2},\dots ,{x}_{n}$, 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=\left[\frac{np}{2}\right]$ 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\le \alpha <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 $\stackrel{^}{f}\left(\omega \right)$ is approximately that of a scaled ${\chi }_{d}^{2}$ 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\left(\omega \right)$ may be constructed as $\left[ml×\stackrel{^}{f}\left(\omega \right),mu×\stackrel{^}{f}\left(\omega \right)\right]$. Alternatively, log $\stackrel{^}{f}\left(\omega \right)$ 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:     $\mathrm{nx}$int64int32nag_int scalar
$n$, the length of the time series.
Constraint: ${\mathbf{nx}}\ge 1$.
2:     $\mathrm{mtx}$int64int32nag_int scalar
If covariances are to be calculated by the function (${\mathbf{ic}}=0$), mtx must specify whether the data are to be initially mean or trend corrected.
${\mathbf{mtx}}=0$
For no correction.
${\mathbf{mtx}}=1$
For mean correction.
${\mathbf{mtx}}=2$
For trend correction.
Constraint: if ${\mathbf{ic}}=0$, $0\le {\mathbf{mtx}}\le 2$
If covariances are supplied (${\mathbf{ic}}\ne 0$), mtx is not used.
3:     $\mathrm{px}$ – double scalar
If covariances are to be calculated by the function (${\mathbf{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 $\left({\mathbf{ic}}\ne 0\right)$, 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\le {\mathbf{px}}\le 1.0$.
4:     $\mathrm{iw}$int64int32nag_int scalar
The choice of lag window.
${\mathbf{iw}}=1$
Rectangular.
${\mathbf{iw}}=2$
Bartlett.
${\mathbf{iw}}=3$
Tukey.
${\mathbf{iw}}=4$
Parzen.
Constraint: $1\le {\mathbf{iw}}\le 4$.
5:     $\mathrm{mw}$int64int32nag_int scalar
$M$, the ‘cut-off’ point of the lag window. Windowed covariances at lag $M$ or greater are zero.
Constraint: $1\le {\mathbf{mw}}\le {\mathbf{nx}}$.
6:     $\mathrm{ic}$int64int32nag_int scalar
Indicates whether covariances are to be calculated in the function or supplied in the call to the function.
${\mathbf{ic}}=0$
Covariances are to be calculated.
${\mathbf{ic}}\ne 0$
Covariances are to be supplied.
7:     $\mathrm{c}\left({\mathbf{nc}}\right)$ – double array
If ${\mathbf{ic}}\ne 0$, c must contain the nc covariances for lags from $0$ to $\left({\mathbf{nc}}-1\right)$, otherwise c need not be set.
8:     $\mathrm{kc}$int64int32nag_int scalar
If ${\mathbf{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 ${2}^{m}$ where $m$ is the smallest integer such that ${2}^{m}\ge {\mathbf{nx}}+{\mathbf{nc}}$, provided $m\le 20$.
If ${\mathbf{ic}}\ne 0$, that is covariances are supplied, kc is not used.
Constraint: ${\mathbf{kc}}\ge {\mathbf{nx}}+{\mathbf{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:     $\mathrm{l}$int64int32nag_int scalar
$L$, the frequency division of the spectral estimates as $\frac{2\pi }{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 ${2}^{m}$ where $m$ is the smallest integer such that ${2}^{m}\ge 2M-1$, provided $m\le 20$.
Constraint: ${\mathbf{l}}\ge 2×{\mathbf{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:   $\mathrm{lg}$int64int32nag_int scalar
Indicates whether unlogged or logged spectral estimates and confidence limits are required.
${\mathbf{lg}}=0$
Unlogged.
${\mathbf{lg}}\ne 0$
Logged.
11:   $\mathrm{xg}\left({\mathbf{nxg}}\right)$ – 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:     $\mathrm{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: ${\mathbf{mw}}\le {\mathbf{nc}}\le {\mathbf{nx}}$.
2:     $\mathrm{nxg}$int64int32nag_int scalar
Default: the dimension of the array xg.
The dimension of the array xg.
Constraints:
• if ${\mathbf{ic}}=0$, ${\mathbf{nxg}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{kc}},{\mathbf{l}}\right)$;
• if ${\mathbf{ic}}\ne 0$, ${\mathbf{nxg}}\ge {\mathbf{l}}$.

### Output Parameters

1:     $\mathrm{c}\left({\mathbf{nc}}\right)$ – double array
If ${\mathbf{ic}}=0$, c will contain the nc calculated covariances.
If ${\mathbf{ic}}\ne 0$, the contents of c will be unchanged.
2:     $\mathrm{xg}\left({\mathbf{nxg}}\right)$ – double array
Contains the ng spectral estimates, $\stackrel{^}{f}\left({\omega }_{\mathit{i}}\right)$, for $\mathit{i}=0,1,\dots ,\left[L/2\right]$ in ${\mathbf{xg}}\left(1\right)$ to ${\mathbf{xg}}\left({\mathbf{ng}}\right)$ respectively (logged if ${\mathbf{lg}}=1$). The elements ${\mathbf{xg}}\left(\mathit{i}\right)$, for $\mathit{i}={\mathbf{ng}}+1,\dots ,{\mathbf{nxg}}$ contain $0.0$.
3:     $\mathrm{ng}$int64int32nag_int scalar
The number of spectral estimates, $\left[L/2\right]+1$, in xg.
4:     $\mathrm{stats}\left(4\right)$ – double array
Four associated statistics. These are the degrees of freedom in ${\mathbf{stats}}\left(1\right)$, the lower and upper $95%$ confidence limit factors in ${\mathbf{stats}}\left(2\right)$ and ${\mathbf{stats}}\left(3\right)$ respectively (logged if ${\mathbf{lg}}=1$), and the bandwidth in ${\mathbf{stats}}\left(4\right)$.
5:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{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.

${\mathbf{ifail}}=1$
 On entry, ${\mathbf{nx}}<1$, or ${\mathbf{mtx}}<0$ and ${\mathbf{ic}}=0$, or ${\mathbf{mtx}}>2$ and ${\mathbf{ic}}=0$, or ${\mathbf{px}}<0.0$, or ${\mathbf{px}}>1.0$, or ${\mathbf{iw}}<1$, or ${\mathbf{iw}}>4$, or ${\mathbf{mw}}<1$, or ${\mathbf{mw}}>{\mathbf{nx}}$, or ${\mathbf{nc}}<{\mathbf{mw}}$, or ${\mathbf{nc}}>{\mathbf{nx}}$, or ${\mathbf{nxg}}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{kc}},{\mathbf{l}}\right)$ and ${\mathbf{ic}}=0$, or ${\mathbf{nxg}}<{\mathbf{l}}$ and ${\mathbf{ic}}\ne 0$.
${\mathbf{ifail}}=2$
 On entry, ${\mathbf{kc}}<{\mathbf{nx}}+{\mathbf{nc}}$, or kc has a prime factor exceeding $19$, or kc has more than $20$ prime factors, counting repetitions.
This error only occurs when ${\mathbf{ic}}=0$.
${\mathbf{ifail}}=3$
 On entry, ${\mathbf{l}}<2×{\mathbf{mw}}-1$, or l has a prime factor exceeding $19$, or l has more than $20$ prime factors, counting repetitions.
W  ${\mathbf{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  ${\mathbf{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.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{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.

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 $n\mathrm{log}\left(n\right)$ (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\pi /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.
```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
```