Mercurial > hg > ltpda

comparison m-toolbox/test/diagnostics/ltpda_arma_time.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 function varargout = ltpda_arma_time(varargin)
2 % LTPDA_ARMA_TIME estimates the ARMA parameters of the transfer function
3 % relating an output to an input
4 %
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 %
7 % DESCRIPTION: LTPDA_ARMA_TIME estimates the ARMA parameters of the
8 % transfer function relating an output to an input time series. The
9 % algorithm uses the IV method to find a set of initial parameters
10 % and then performs an iterative search based on the Optimization
11 % Toolbox.
12 %
13 % CALL: b = ltpda_arma_time(ax,ay,pl)
14 %
15 % INPUTS: ax - analysis object containig the input time series
16 % ay - analysis object containig the output time series
17 % pl - parameters list
18 %
19 % OUTPUTS: b - cdata type analysis object containing the ARMA
20 % parameters and the residual of the fit
21 %
22 % PARAMETERS:
23 % MaxIter - Maximum number of iterations to be performed by the
24 % iterative search (default 4e2)
25 % MaxFunEvals - Maximum number of function evaluations (default 1e2)
26 % TolX - Tolerance on the estimated parameter (default 1e-6)
27 % TolFun - Tolerance on the evaluated function (default 1e-6)
28 % UpBound - Array of parameters upper bounds for the iterative search
29 % (default is 1 for each parameter)
30 % LowBound - Array of parameters lower bounds for the iterative search
31 % (default is 1 for each parameter)
32 % ARMANum - MA order of the ARMA filter (default 2)
33 % ARMADen - AR order of the ARMA filter (default 1)
34 %
35 % VERSION: $Id: ltpda_arma_time.m,v 1.1 2008/03/04 12:55:05 miquel Exp $
36 %
37 % HISTORY: 15-02-2008 M Nofrarias
38 % Creation
39 %
40 %
41 % TODO: - Add parameters errors
42 % - Pass from univariate to multivariate
43 %
44 %
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46
47
48
49 %% Standard history variables
50
51 ALGONAME = mfilename;
52 VERSION = '$Id: ltpda_arma_time.m,v 1.1 2008/03/04 12:55:05 miquel Exp $';
53 CATEGORY = 'SIGPROC';
54
55
56 %% Check if this is a call for parameters
57 if nargin == 1 && ischar(varargin{1})
58 in = char(varargin{1});
59 if strcmpi(in, 'Params')
60 varargout{1} = getDefaultPL();
61 return
62 elseif strcmpi(in, 'Version')
63 varargout{1} = VERSION;
64 return
65 elseif strcmpi(in, 'Category')
66 varargout{1} = CATEGORY;
67 return
68 end
69 end
70
71 %% Capture input variables names
72 invars = {};
73 as = [];
74 ps = [];
75 for j=1:nargin
76 invars = [invars cellstr(inputname(j))];
77 if isa(varargin{j}, 'ao')
78 as = [as varargin{j}];
79 end
80 if isa(varargin{j}, 'plist')
81 ps = [ps varargin{j}];
82 end
83 end
84
85
86 %% Check plist
87 if isempty(ps)
88 pl = getDefaultPL();
89 else
90 pl = combine(ps, getDefaultPL);
91 end
92
93
94 %% Initialise variables
95 x=as(1).data.y;
96 y=as(2).data.y;
97 p = find(pl, 'ARMANum');
98 q = find(pl, 'ARMADen');
99 disp(sprintf('! ARMA orders p = %d, q = %d', p,q))
100
101 %% First estimate through Instrument Variables (IV) Method
102
103 % Least squares
104 for i=1:p
105 obs1(:,i) = [zeros(i,1); x(1:length(x)-i)];
106 end
107 for i=1:q
108 obs2(:,i) = [zeros(i,1); -y(1:length(y)-i)];
109 end
110
111 % LS estimation
112 obs = [obs1 obs2];
113 par0=(obs'*obs)\obs'*y
114
115 % filter observations with LS parameters
116 obs1_f = filter([par0(1:p)],[1 par0(p+1:p+q)],obs1);
117
118 % build instruments
119 inst = [obs1 obs1_f];
120
121 % IV estimation
122 par0=(inst'*obs)\inst'*y
123
124
125 %% Iterative search using the Optimization Toolbox
126
127
128 % Options for the lsqcurvefit function
129 opt = optimset(...
130 'MaxIter',find(pl, 'MaxIter'),...
131 'TolX',find(pl, 'TolXm'),...
132 'TolFun',find(pl, 'TolFun'),...
133 'MaxFunEvals',find(pl, 'MaxFunEvals'),...
134 'Display','off');
135
136 % Upper and Lower Bounds for the parameters
137 ub=find(pl,'UpBound');
138 lb=find(pl,'LowBound');
139 if isempty(ub)
140 ub=ones(p+q,1);
141 pl = append(pl, param('UpBound', ub));
142 disp('! Parameters Upper Bound set to 1')
143 end
144 if isempty(lb)
145 lb=-ones(p+q,1);
146 pl = append(pl, param('LowBound', lb));
147 disp('! Parameters Lower Bound set to -1')
148 end
149
150 % Call to the Optimization Toolbox function
151 [par,res]=lsqcurvefit(@(par,xdata)filter([par(1:p)],[1 par(p+1:p+q)],xdata),...
152 par0,x,y,lb,ub,opt);
153
154
155 %% Build output ao
156
157 % New output history
158 h = history(ALGONAME, VERSION, pl,[as(1).hist as(2).hist]);
159 h = set(h,'invars', invars);
160
161 % Make output analysis object
162 b = ao([par; res]);
163
164 % Name, mfilename, description for this object
165 b = setnh(b,...
166 'name', ['ARMA param. Input: ' sprintf('%s ', char(invars{1})) ' Output: ' char(invars{2})],...
167 'mfilename', ALGONAME, ...
168 'description', find(pl,'description'),...
169 'hist', h);
170
171 varargout{1} = b;
172
173
174 %% --- FUNCTIONS ---
175
176
177 % Get default parameters
178 function plo = getDefaultPL();
179 disp('* creating default plist...');
180 plo = plist();
181 plo = append(plo, param('MaxIter',4e2));
182 plo = append(plo, param('MaxFunEvals',1e2));
183 plo = append(plo, param('TolX',1e-6));
184 plo = append(plo, param('TolFun',1e-6));
185 plo = append(plo, param('ARMANum',2));
186 plo = append(plo, param('ARMADen',1));
187 disp('* done.');