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
(ω) = 1/(2π)
 ( M − 1 ) C0 + 2 ∑ wkCkcos(ωk) k = 1
,
$f^(ω)=12π (C0+2∑k=1 M-1wkCkcos(ωk)) ,$
where M$M$ is the window width, and is calculated for frequency values
 ωi = (2πi)/L,  i = 0,1, … ,[L / 2], $ωi=2πiL, i=0,1,…,[L/2],$
where []$\left[\right]$ denotes the integer part.
The autocovariances Ck${C}_{k}$ may be supplied by you, or constructed from a time series x1,x2,,xn${x}_{1},{x}_{2},\dots ,{x}_{n}$, as
 n − k Ck = 1/n ∑ xtxt + k, t = 1
$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:
 (1/2)(1 − cos(π(t − (1/2)) / T)), 1 ≤ t ≤ T (1/2)(1 − cos(π(n − t + (1/2)) / T)), n + 1 − T ≤ t ≤ n 1, otherwise,
$12(1-cos(π (t-12)/T)), 1≤t≤T 12(1-cos(π (n-t+12)/T)), n+ 1-T≤t≤n 1, otherwise,$
where T = [(np)/2] $T=\left[\frac{np}{2}\right]$ and p$p$ is the tapering proportion.
The smoothing window is defined by
 wk = W (k/M) ,  k ≤ M − 1, $wk=W (kM) , k≤M-1,$
which for the various windows is defined over 0α < 1$0\le \alpha <1$ by
rectangular:
 W(α) = 1 $W(α)=1$
Bartlett:
 W(α) = 1 − α $W(α )= 1-α$
Tukey:
 W(α) = (1/2)(1 + cos(πα)) $W(α)=12(1+cos(πα))$
Parzen:
 W(α) = 1 − 6α2 + 6α3, 0 ≤ α ≤ (1/2) W(α) = 2(1 − α)3, (1/2) < α < 1.
$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 χd2${\chi }_{d}^{2}$ variate, whose degrees of freedom d$d$ is provided by the function, together with multiplying limits mu$mu$, ml$ml$ from which approximate 95%$95%$ confidence intervals for the true spectrum f(ω)$f\left(\omega \right)$ may be constructed as [ ml × (ω) , mu × (ω) ] $\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$b$ of the corresponding smoothing window in the frequency domain is also provided. Spectrum estimates separated by (angular) frequencies much greater than b$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$n$, the length of the time series.
Constraint: nx1${\mathbf{nx}}\ge 1$.
2:     mtx – int64int32nag_int scalar
If covariances are to be calculated by the function (ic = 0${\mathbf{ic}}=0$), mtx must specify whether the data are to be initially mean or trend corrected.
mtx = 0${\mathbf{mtx}}=0$
For no correction.
mtx = 1${\mathbf{mtx}}=1$
For mean correction.
mtx = 2${\mathbf{mtx}}=2$
For trend correction.
Constraint: if ic = 0${\mathbf{ic}}=0$, 0mtx2$0\le {\mathbf{mtx}}\le 2$
If covariances are supplied (ic0${\mathbf{ic}}\ne 0$), mtx is not used.
3:     px – double scalar
If covariances are to be calculated by the function (ic = 0${\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 (ic0)$\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$0.0$ implies no tapering.
Constraint: 0.0px1.0$0.0\le {\mathbf{px}}\le 1.0$.
4:     iw – int64int32nag_int scalar
The choice of lag window.
iw = 1${\mathbf{iw}}=1$
Rectangular.
iw = 2${\mathbf{iw}}=2$
Bartlett.
iw = 3${\mathbf{iw}}=3$
Tukey.
iw = 4${\mathbf{iw}}=4$
Parzen.
Constraint: 1iw4$1\le {\mathbf{iw}}\le 4$.
5:     mw – int64int32nag_int scalar
M$M$, the ‘cut-off’ point of the lag window. Windowed covariances at lag M$M$ or greater are zero.
Constraint: 1mwnx$1\le {\mathbf{mw}}\le {\mathbf{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${\mathbf{ic}}=0$
Covariances are to be calculated.
ic0${\mathbf{ic}}\ne 0$
Covariances are to be supplied.
7:     c(nc) – double array
nc, the dimension of the array, must satisfy the constraint mwncnx${\mathbf{mw}}\le {\mathbf{nc}}\le {\mathbf{nx}}$.
If ic0${\mathbf{ic}}\ne 0$, c must contain the nc covariances for lags from 0$0$ to (nc1)$\left({\mathbf{nc}}-1\right)$, otherwise c need not be set.
8:     kc – int64int32nag_int scalar
If ic = 0${\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 2m${2}^{m}$ where m$m$ is the smallest integer such that 2mnx + nc${2}^{m}\ge {\mathbf{nx}}+{\mathbf{nc}}$, provided m20$m\le 20$.
If ic0${\mathbf{ic}}\ne 0$, that is covariances are supplied, kc is not used.
Constraint: kcnx + nc${\mathbf{kc}}\ge {\mathbf{nx}}+{\mathbf{nc}}$. The largest prime factor of kc must not exceed 19$19$, and the total number of prime factors of kc, counting repetitions, must not exceed 20$20$. These two restrictions are imposed by the internal FFT algorithm used.
9:     l – int64int32nag_int scalar
L$L$, the frequency division of the spectral estimates as (2π)/L $\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 2m${2}^{m}$ where m$m$ is the smallest integer such that 2m2M1${2}^{m}\ge 2M-1$, provided m20$m\le 20$.
Constraint: l2 × mw1${\mathbf{l}}\ge 2×{\mathbf{mw}}-1$. The largest prime factor of l must not exceed 19$19$, and the total number of prime factors of l, counting repetitions, must not exceed 20$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${\mathbf{lg}}=0$
Unlogged.
lg0${\mathbf{lg}}\ne 0$
Logged.
11:   xg(nxg) – double array
nxg, the dimension of the array, must satisfy the constraint
• if ic = 0${\mathbf{ic}}=0$, nxgmax (kc,l)${\mathbf{nxg}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{kc}},{\mathbf{l}}\right)$;
• if ic0${\mathbf{ic}}\ne 0$, nxgl${\mathbf{nxg}}\ge {\mathbf{l}}$.
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: mwncnx${\mathbf{mw}}\le {\mathbf{nc}}\le {\mathbf{nx}}$.
2:     nxg – int64int32nag_int scalar
Default: The dimension of the array xg.
The dimension of the array xg as declared in the (sub)program from which nag_tsa_uni_spectrum_lag (g13ca) is called.
Constraints:
• if ic = 0${\mathbf{ic}}=0$, nxgmax (kc,l)${\mathbf{nxg}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{kc}},{\mathbf{l}}\right)$;
• if ic0${\mathbf{ic}}\ne 0$, nxgl${\mathbf{nxg}}\ge {\mathbf{l}}$.

None.

### Output Parameters

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

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

## 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$n$ is approximately proportional to nlog(n)$n\mathrm{log}\left(n\right)$ (but see Section [Further Comments] in (c06pa) for further details).

## Example

```function nag_tsa_uni_spectrum_lag_example
nx = int64(256);
mtx = int64(1);
px = 0.1;
iw = int64(4);
mw = int64(100);
ic = int64(0);
c = zeros(100, 1);
kc = int64(360);
l = int64(200);
lg = int64(0);
xg = zeros(360,1);
xg(1:256) = [5;
11;
16;
23;
36;
58;
29;
20;
10;
8;
3;
0;
0;
2;
11;
27;
47;
63;
60;
39;
28;
26;
22;
11;
21;
40;
78;
122;
103;
73;
47;
35;
11;
5;
16;
34;
70;
81;
111;
101;
73;
40;
20;
16;
5;
11;
22;
40;
60;
80.9;
83.4;
47.7;
47.8;
30.7;
12.2;
9.6;
10.2;
32.4;
47.6;
54;
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;
19.8;
92.5;
154.4;
125.9;
84.8;
68.1;
38.5;
22.8;
10.2;
24.1;
82.9;
132;
130.9;
118.1;
89.9;
66.6;
60;
46.9;
41;
21.3;
16;
6.4;
4.1;
6.8;
14.5;
34;
45;
43.1;
47.5;
42.2;
28.1;
10.1;
8.1;
2.5;
0;
1.4;
5;
12.2;
13.9;
35.4;
45.8;
41.1;
30.1;
23.9;
15.6;
6.6;
4;
1.8;
8.5;
16.6;
36.3;
49.6;
64.2;
67;
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;
40.1;
61.5;
98.5;
124.7;
96.3;
66.6;
64.5;
54.1;
39;
20.6;
6.7;
4.3;
22.7;
54.8;
93.8;
95.8;
77.2;
59.1;
44;
47;
30.5;
16.3;
7.3;
37.6;
74;
139;
111.2;
101.6;
66.2;
44.7;
17;
11.3;
12.4;
3.4;
6;
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;
85.1;
78;
64;
41.8;
26.2;
26.7;
12.1;
9.5;
2.7;
5;
24.4;
42;
63.5;
53.8;
62;
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;
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];
[cOut, xgOut, ng, stats, ifail] = ...
nag_tsa_uni_spectrum_lag(nx, mtx, px, iw, mw, ic, c, kc, l, lg, xg);
ng, stats, ifail
```
```

ng =

101

stats =

9.0000
0.4731
3.3329
0.1165

ifail =

0

```
```function g13ca_example
nx = int64(256);
mtx = int64(1);
px = 0.1;
iw = int64(4);
mw = int64(100);
ic = int64(0);
c = zeros(100, 1);
kc = int64(360);
l = int64(200);
lg = int64(0);
xg = zeros(360,1);
xg(1:256) = [5;
11;
16;
23;
36;
58;
29;
20;
10;
8;
3;
0;
0;
2;
11;
27;
47;
63;
60;
39;
28;
26;
22;
11;
21;
40;
78;
122;
103;
73;
47;
35;
11;
5;
16;
34;
70;
81;
111;
101;
73;
40;
20;
16;
5;
11;
22;
40;
60;
80.9;
83.4;
47.7;
47.8;
30.7;
12.2;
9.6;
10.2;
32.4;
47.6;
54;
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;
19.8;
92.5;
154.4;
125.9;
84.8;
68.1;
38.5;
22.8;
10.2;
24.1;
82.9;
132;
130.9;
118.1;
89.9;
66.6;
60;
46.9;
41;
21.3;
16;
6.4;
4.1;
6.8;
14.5;
34;
45;
43.1;
47.5;
42.2;
28.1;
10.1;
8.1;
2.5;
0;
1.4;
5;
12.2;
13.9;
35.4;
45.8;
41.1;
30.1;
23.9;
15.6;
6.6;
4;
1.8;
8.5;
16.6;
36.3;
49.6;
64.2;
67;
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;
40.1;
61.5;
98.5;
124.7;
96.3;
66.6;
64.5;
54.1;
39;
20.6;
6.7;
4.3;
22.7;
54.8;
93.8;
95.8;
77.2;
59.1;
44;
47;
30.5;
16.3;
7.3;
37.6;
74;
139;
111.2;
101.6;
66.2;
44.7;
17;
11.3;
12.4;
3.4;
6;
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;
85.1;
78;
64;
41.8;
26.2;
26.7;
12.1;
9.5;
2.7;
5;
24.4;
42;
63.5;
53.8;
62;
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;
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];
[cOut, xgOut, ng, stats, ifail] = g13ca(nx, mtx, px, iw, mw, ic, c, kc, l, lg, xg);
ng, stats, ifail
```
```

ng =

101

stats =

9.0000
0.4731
3.3329
0.1165

ifail =

0

```