LTPDA_ARMA_TIME estimates the ARMA parameters of the transfer function relating an output to an input %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION: LTPDA_ARMA_TIME estimates the ARMA parameters of the transfer function relating an output to an input time series. The algorithm uses the IV method to find a set of initial parameters and then performs an iterative search based on the Optimization Toolbox. CALL: b = ltpda_arma_time(ax,ay,pl) INPUTS: ax - analysis object containig the input time series ay - analysis object containig the output time series pl - parameters list OUTPUTS: b - cdata type analysis object containing the ARMA parameters and the residual of the fit PARAMETERS: MaxIter - Maximum number of iterations to be performed by the iterative search (default 4e2) MaxFunEvals - Maximum number of function evaluations (default 1e2) TolX - Tolerance on the estimated parameter (default 1e-6) TolFun - Tolerance on the evaluated function (default 1e-6) UpBound - Array of parameters upper bounds for the iterative search (default is 1 for each parameter) LowBound - Array of parameters lower bounds for the iterative search (default is 1 for each parameter) ARMANum - MA order of the ARMA filter (default 2) ARMADen - AR order of the ARMA filter (default 1) VERSION: $Id: ltpda_arma_time.html,v 1.1 2008/03/01 12:29:33 hewitson Exp $ HISTORY: 15-02-2008 M Nofrarias Creation TODO: - Add parameters errors - Pass from univariate to multivariate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 function varargout = ltpda_arma_time(varargin) 0002 % LTPDA_ARMA_TIME estimates the ARMA parameters of the transfer function 0003 % relating an output to an input 0004 % 0005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0006 % 0007 % DESCRIPTION: LTPDA_ARMA_TIME estimates the ARMA parameters of the 0008 % transfer function relating an output to an input time series. The 0009 % algorithm uses the IV method to find a set of initial parameters 0010 % and then performs an iterative search based on the Optimization 0011 % Toolbox. 0012 % 0013 % CALL: b = ltpda_arma_time(ax,ay,pl) 0014 % 0015 % INPUTS: ax - analysis object containig the input time series 0016 % ay - analysis object containig the output time series 0017 % pl - parameters list 0018 % 0019 % OUTPUTS: b - cdata type analysis object containing the ARMA 0020 % parameters and the residual of the fit 0021 % 0022 % PARAMETERS: 0023 % MaxIter - Maximum number of iterations to be performed by the 0024 % iterative search (default 4e2) 0025 % MaxFunEvals - Maximum number of function evaluations (default 1e2) 0026 % TolX - Tolerance on the estimated parameter (default 1e-6) 0027 % TolFun - Tolerance on the evaluated function (default 1e-6) 0028 % UpBound - Array of parameters upper bounds for the iterative search 0029 % (default is 1 for each parameter) 0030 % LowBound - Array of parameters lower bounds for the iterative search 0031 % (default is 1 for each parameter) 0032 % ARMANum - MA order of the ARMA filter (default 2) 0033 % ARMADen - AR order of the ARMA filter (default 1) 0034 % 0035 % VERSION: $Id: ltpda_arma_time.html,v 1.1 2008/03/01 12:29:33 hewitson Exp $ 0036 % 0037 % HISTORY: 15-02-2008 M Nofrarias 0038 % Creation 0039 % 0040 % 0041 % TODO: - Add parameters errors 0042 % - Pass from univariate to multivariate 0043 % 0044 % 0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0046 0047 0048 0049 %% Standard history variables 0050 0051 ALGONAME = mfilename; 0052 VERSION = '$Id: ltpda_arma_time.html,v 1.1 2008/03/01 12:29:33 hewitson Exp $'; 0053 CATEGORY = 'SIGPROC'; 0054 0055 0056 %% Check if this is a call for parameters 0057 if nargin == 1 && ischar(varargin{1}) 0058 in = char(varargin{1}); 0059 if strcmpi(in, 'Params') 0060 varargout{1} = getDefaultPL(); 0061 return 0062 elseif strcmpi(in, 'Version') 0063 varargout{1} = VERSION; 0064 return 0065 elseif strcmpi(in, 'Category') 0066 varargout{1} = CATEGORY; 0067 return 0068 end 0069 end 0070 0071 %% Capture input variables names 0072 invars = {}; 0073 as = []; 0074 ps = []; 0075 for j=1:nargin 0076 invars = [invars cellstr(inputname(j))]; 0077 if isa(varargin{j}, 'ao') 0078 as = [as varargin{j}]; 0079 end 0080 if isa(varargin{j}, 'plist') 0081 ps = [ps varargin{j}]; 0082 end 0083 end 0084 0085 0086 %% Check plist 0087 if isempty(ps) 0088 pl = getDefaultPL(); 0089 else 0090 pl = combine(ps, getDefaultPL); 0091 end 0092 0093 0094 %% Initialise variables 0095 x=as(1).data.y; 0096 y=as(2).data.y; 0097 p = find(pl, 'ARMANum'); 0098 q = find(pl, 'ARMADen'); 0099 disp(sprintf('! ARMA orders p = %d, q = %d', p,q)) 0100 0101 %% First estimate through Instrument Variables (IV) Method 0102 0103 for i=1:p 0104 inst1(:,i) = [zeros(i,1); x(1:length(x)-i)]; 0105 end 0106 for i=1:q 0107 inst2(:,i) = [zeros(i,1); -y(1:length(y)-i)]; 0108 end 0109 0110 inst = [inst1 inst2]; 0111 par0=(inst'*inst)\inst'*y; 0112 0113 %% Iterative search using the Optimization Toolbox 0114 0115 0116 % Options for the lsqcurvefit function 0117 opt = optimset(... 0118 'MaxIter',find(pl, 'MaxIter'),... 0119 'TolX',find(pl, 'TolXm'),... 0120 'TolFun',find(pl, 'TolFun'),... 0121 'MaxFunEvals',find(pl, 'MaxFunEvals'),... 0122 'Display','off'); 0123 0124 % Upper and Lower Bounds for the parameters 0125 ub=find(pl,'UpBound'); 0126 lb=find(pl,'LowBound'); 0127 if isempty(ub) 0128 ub=ones(p+q,1); 0129 pl = append(pl, param('UpBound', ub)); 0130 disp('! Parameters Upper Bound set to 1') 0131 end 0132 if isempty(lb) 0133 lb=-ones(p+q,1); 0134 pl = append(pl, param('LowBound', lb)); 0135 disp('! Parameters Lower Bound set to -1') 0136 end 0137 0138 % Call to the Optimization Toolbox function 0139 [par,res]=lsqcurvefit(@(par,xdata)filter([par(1:p)],[1 par(p+1:p+q)],xdata),... 0140 par0,x,y,lb,ub,opt); 0141 0142 0143 %% Build output ao 0144 0145 % New output history 0146 h = history(ALGONAME, VERSION, pl,[as(1).hist as(2).hist]); 0147 h = set(h,'invars', invars); 0148 0149 % Make output analysis object 0150 b = ao([par; res]); 0151 0152 % Name, mfilename, description for this object 0153 b = setnh(b,... 0154 'name', ['ARMA param. Input: ' sprintf('%s ', char(invars{1})) ' Output: ' char(invars{2})],... 0155 'mfilename', ALGONAME, ... 0156 'description', find(pl,'description'),... 0157 'hist', h); 0158 0159 varargout{1} = b; 0160 0161 0162 %% --- FUNCTIONS --- 0163 0164 0165 % Get default parameters 0166 function plo = getDefaultPL(); 0167 disp('* creating default plist...'); 0168 plo = plist(); 0169 plo = append(plo, param('MaxIter',4e2)); 0170 plo = append(plo, param('MaxFunEvals',1e2)); 0171 plo = append(plo, param('TolX',1e-6)); 0172 plo = append(plo, param('TolFun',1e-6)); 0173 plo = append(plo, param('ARMANum',2)); 0174 plo = append(plo, param('ARMADen',1)); 0175 disp('* done.');