Mercurial > octave
annotate scripts/signal/arch_fit.m @ 31225:3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
| author | Arun Giridhar <arungiridhar@gmail.com> |
|---|---|
| date | Sun, 11 Sep 2022 13:53:38 -0400 |
| parents | 796f54d4ddbf |
| children | 597f3ee61a48 |
| rev | line source |
|---|---|
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 ######################################################################## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 ## |
|
30564
796f54d4ddbf
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
29359
diff
changeset
|
3 ## Copyright (C) 1995-2022 The Octave Project Developers |
|
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
4 ## |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 ## See the file COPYRIGHT.md in the top-level directory of this |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 ## distribution or <https://octave.org/copyright/>. |
| 3426 | 7 ## |
| 3922 | 8 ## This file is part of Octave. |
| 9 ## | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
10 ## Octave is free software: you can redistribute it and/or modify it |
| 3922 | 11 ## under the terms of the GNU General Public License as published by |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
12 ## the Free Software Foundation, either version 3 of the License, or |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
13 ## (at your option) any later version. |
| 3426 | 14 ## |
| 3922 | 15 ## Octave is distributed in the hope that it will be useful, but |
| 3191 | 16 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
18 ## GNU General Public License for more details. |
| 3426 | 19 ## |
| 3191 | 20 ## You should have received a copy of the GNU General Public License |
| 7016 | 21 ## along with Octave; see the file COPYING. If not, see |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
22 ## <https://www.gnu.org/licenses/>. |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 ## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 ######################################################################## |
| 3191 | 25 |
| 3449 | 26 ## -*- texinfo -*- |
|
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
27 ## @deftypefn {} {[@var{a}, @var{b}] =} arch_fit (@var{y}, @var{x}, @var{p}, @var{iter}, @var{gamma}, @var{a0}, @var{b0}) |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
28 ## Fit an ARCH regression model to the time series @var{y} using the scoring |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
29 ## algorithm in @nospell{Engle's} original ARCH paper. |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
30 ## |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
31 ## The model is |
| 3191 | 32 ## |
| 3449 | 33 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
34 ## @group |
|
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
35 ## y(t) = b(1) * x(t,1) + @dots{} + b(k) * x(t,k) + e(t), |
|
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
36 ## h(t) = a(1) + a(2) * e(t-1)^2 + @dots{} + a(p+1) * e(t-p)^2 |
|
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
37 ## @end group |
| 3449 | 38 ## @end example |
| 3191 | 39 ## |
| 3449 | 40 ## @noindent |
| 3499 | 41 ## in which @math{e(t)} is @math{N(0, h(t))}, given a time-series vector |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
42 ## @var{y} up to time @math{t-1} and a matrix of (ordinary) regressors @var{x} |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
43 ## up to @math{t}. The order of the regression of the residual variance is |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
44 ## specified by @var{p}. |
| 3191 | 45 ## |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
46 ## If invoked as @code{arch_fit (@var{y}, @var{k}, @var{p})} with a positive |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
47 ## integer @var{k}, fit an ARCH(@var{k}, @var{p}) process, i.e., do the above |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
48 ## with the @math{t}-th row of @var{x} given by |
| 3191 | 49 ## |
| 3449 | 50 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
51 ## [1, y(t-1), @dots{}, y(t-k)] |
| 3449 | 52 ## @end example |
| 53 ## | |
| 54 ## Optionally, one can specify the number of iterations @var{iter}, the | |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
55 ## updating factor @var{gamma}, and initial values @math{a0} and @math{b0} |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
56 ## for the scoring algorithm. |
| 3449 | 57 ## @end deftypefn |
| 3191 | 58 |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
59 function [a, b] = arch_fit (y, x, p, iter, gamma, a0, b0) |
| 3191 | 60 |
|
28789
28de41192f3c
Eliminate unneeded verification of nargin, nargout in m-files.
Rik <rik@octave.org>
parents:
27984
diff
changeset
|
61 if (nargin < 3 || nargin == 6) |
| 6046 | 62 print_usage (); |
| 3191 | 63 endif |
| 3426 | 64 |
| 4030 | 65 if (! (isvector (y))) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11469
diff
changeset
|
66 error ("arch_fit: Y must be a vector"); |
| 3191 | 67 endif |
| 68 | |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
69 T = length (y); |
|
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
70 y = reshape (y, T, 1); |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
71 [rx, cx] = size (x); |
| 3191 | 72 if ((rx == 1) && (cx == 1)) |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
73 x = autoreg_matrix (y, x); |
| 3457 | 74 elseif (! (rx == T)) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11469
diff
changeset
|
75 error ("arch_fit: either rows (X) == length (Y), or X is a scalar"); |
| 3191 | 76 endif |
| 77 | |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
78 [T, k] = size (x); |
| 3426 | 79 |
| 3191 | 80 if (nargin == 7) |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
81 a = a0; |
|
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
82 b = b0; |
|
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
83 e = y - x * b; |
| 3191 | 84 else |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
85 [b, v_b, e] = ols (y, x); |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
86 a = [v_b, (zeros (1, p))]'; |
| 3191 | 87 if (nargin < 5) |
| 88 gamma = 0.1; | |
| 89 if (nargin < 4) | |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
90 iter = 50; |
| 3191 | 91 endif |
| 92 endif | |
| 93 endif | |
| 3426 | 94 |
| 3191 | 95 esq = e.^2; |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
96 Z = autoreg_matrix (esq, p); |
| 3191 | 97 |
|
26268
6dd232798997
maint: Remove useless ';' from end of for, if, while, etc. statements.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
98 for i = 1 : iter |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
99 h = Z * a; |
|
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
100 tmp = esq ./ h.^2 - 1 ./ h; |
|
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
101 s = 1 ./ h(1:T-p); |
|
26268
6dd232798997
maint: Remove useless ';' from end of for, if, while, etc. statements.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
102 for j = 1 : p |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
103 s -= a(j+1) * tmp(j+1:T-p+j); |
| 3191 | 104 endfor |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
105 r = 1 ./ h(1:T-p); |
|
26268
6dd232798997
maint: Remove useless ';' from end of for, if, while, etc. statements.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
106 for j = 1:p |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
107 r += 2 * h(j+1:T-p+j).^2 .* esq(1:T-p); |
| 3191 | 108 endfor |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
109 r = sqrt (r); |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
110 X_tilde = x(1:T-p, :) .* (r * ones (1,k)); |
| 3191 | 111 e_tilde = e(1:T-p) .*s ./ r; |
| 112 delta_b = inv (X_tilde' * X_tilde) * X_tilde' * e_tilde; | |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
113 b += gamma * delta_b; |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
114 e = y - x * b; |
| 3191 | 115 esq = e .^ 2; |
| 116 Z = autoreg_matrix (esq, p); | |
| 117 h = Z * a; | |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
118 f = esq ./ h - ones (T,1); |
| 3191 | 119 Z_tilde = Z ./ (h * ones (1, p+1)); |
| 120 delta_a = inv (Z_tilde' * Z_tilde) * Z_tilde' * f; | |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
121 a += gamma * delta_a; |
| 3191 | 122 endfor |
| 3426 | 123 |
| 3191 | 124 endfunction |