Added arima class with estimate method.

Use X13 Census for estimation.
parent 8cbcc968
classdef arima<handle % --*-- Unitary tests --*--
% Class for ARIMA models (interface to X13).
% Copyright (C) 2017 Dynare Team
%
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare dates submodule is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
properties
y = []; % dseries object with a single variable.
p = []; % Number of lags on the autoregressive part.
d = []; % Order of differenciation.
q = []; % Number of lags on the moving average part.
P = []; % Seasonal component autoregressive part number of lags.
D = []; % Seasonal component order of differenciation.
Q = []; % Seasonal component moving average part number of lags.
S = []; % Seasonal component period.
ar = []; % Autoregressive parameters.
ma = []; % Moving arverage parameters.
AR = []; % Seasonal component autoregressive parameters.
MA = []; % Seasonal component moving average parameters.
ar_std = []; % Autoregressive parameters estimator std.
ma_std = []; % Moving arverage parameters estimator std.
AR_std = []; % Seasonal component autoregressive parameters estimator std.
MA_std = []; % Seasonal component moving average parameters estimator std.
sigma = []; % Standard deviation of the schock.
estimation = []; % Structure gathering informations related to the estimation of the model.
end
methods
function o = arima(y, model)
% Constructor for the arma class.
%
% INPUTS
% - y [dseries] Data.
% - model [string] Specification of the model of the form '(p, d, q)(P,D,Q)S', See X-13ARIMA-SEATS manual.
%
% OUPUTS
% - o [arima] ARIMA model object.
nargin
switch nargin
case 0
% Return empty object.
o.y = []; % dseries object with a single variable.
o.p = []; % Number of lags on the autoregressive part.
o.d = []; % Order of differenciation.
o.q = []; % Number of lags on the moving average part.
o.P = []; % Seasonal component autoregressive part number of lags.
o.D = []; % Seasonal component order of differenciation.
o.Q = []; % Seasonal component moving average part number of lags.
o.S = []; % Seasonal component period.
o.ar = []; % Autoregressive parameters.
o.ma = []; % Moving arverage parameters.
o.AR = []; % Seasonal component autoregressive parameters.
o.MA = []; % Seasonal component moving average parameters.
o.ar_std = []; % Autoregressive parameters estimator std.
o.ma_std = []; % Moving arverage parameters estimator std.
o.AR_std = []; % Seasonal component autoregressive parameters estimator std.
o.MA_std = []; % Seasonal component moving average parameters estimator std.
o.sigma = []; % Standard deviation of the schock.
o.estimation = []; % Structure gathering informations related to the estimation of the model.
return
case 2
if ~isdseries(y)
error('arima::WrongInputArguments', 'First input argument must be a dseries object!')
end
if ~ischar(model)
error('arima::WrongInputArguments', 'Second input argument must be a string!')
end
o.y = y;
o.estimation = struct();
% Read the description of the ARIMA model.
model = regexprep(model,'(\s)',''); % Removes all spaces.
description = regexp(model, '(\([0-9],[0-9],[0-9]\)([0-9]+)?)', 'tokens');
% First cell is the non seasonal component.
ns_arima = description{1}{1};
ns_arima_specification = regexp(ns_arima, '([0-9]*)', 'tokens');
o.p = str2num(ns_arima_specification{1}{1});
o.d = str2num(ns_arima_specification{2}{1});
o.q = str2num(ns_arima_specification{3}{1});
% Set default values for parameters and estimate std.
if o.p
o.ar = zeros(o.p, 1);
o.ar_std = zeros(o.p, 1);
else
o.ar = [];
o.ar_std = [];
end
if o.q
o.ma = zeros(o.q, 1);
o.ma_std = zeros(o.q, 1);
else
o.ma = [];
o.ma_std = [];
end
% Following cells are the ARIMA models on the seasonal components.
number_of_seasonal_components = length(description)-1;
if number_of_seasonal_components
o.P = zeros(number_of_seasonal_components, 1);
o.D = zeros(number_of_seasonal_components, 1);
o.Q = zeros(number_of_seasonal_components, 1);
o.S = zeros(number_of_seasonal_components, 1);
for i=1:number_of_seasonal_components
s_arima = description{i+1}{1};
s_arima_specification = regexp(s_arima, '([0-9]*)', 'tokens');
if isequal(i, 1) && isequal(length(s_arima_specification), 3)
o.S(i) = y.init.freq;
else
o.S(i) = str2num(s_arima_specification{4}{1});
end
o.P(i) = str2num(s_arima_specification{1}{1});
o.D(i) = str2num(s_arima_specification{2}{1});
o.Q(i) = str2num(s_arima_specification{3}{1});
end
end
% Set default values for parameters and estimate std.
o.AR = {};
o.MA = {};
o.AR_std = {};
o.MA_std = {};
for i=1:number_of_seasonal_components
if o.P(i)
o.AR{i} = zeros(o.P(i), 1);
o.AR_std{i} = zeros(o.P(i), 1);
else
o.AR{i} = [];
o.AR_std{i} = [];
end
if o.Q(i)
o.MA{i} = zeros(o.Q(i), 1);
o.MA_std{i} = zeros(o.Q(i), 1);
else
o.MA{i} = [];
o.MA_std{i} = [];
end
end
% Set default value for the size of the innovation.
o.sigma = std(y.data);
otherwise
error('arima::WrongInputArguments', 'Two input argument are mandatory!')
end
end % arima
end % methods
end % classdef
%@test:1
%$ try
%$ data = dseries(cumsum(randn(200, 1)), '1990Q1', 'Output');
%$ arma = arima(data, '(1,1,1)(0,0,1)');
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ T = all(t);
%@eof:1
function estimate(o, varargin)
% Estimate method for arima class.
% Copyright (C) 2017 Dynare Team
%
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare dates submodule is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Set default values for options.
options.tol = 1e-5;
options.exact = 'arma';
options.maxiter = 1500;
options.outofsample = false;
options.print='none';
options.save = 'estimates lkstats residuals';
% Update options.
if length(varargin)
if mod(length(varargin), 2)
error('arima:estimate:WrongInputArguments', 'Number of input should be even (option name, option value)!')
else
for option_id = 1:2:length(varargin)
if ismember(varargin{option_id}, {'tol', 'exact', 'maxiter', 'outofsample', 'print', 'save'})
eval(sprintf('options.%s = %s;', varargin{option_id}, varargin{option_id+1}))
else
warning('arima:estimate: Unknown option %s!', varargin{option_id})
end
end
end
end
model = sprintf('(%i,%i,%i)', o.p, o.d, o.q);
for i = 1:length(o.P)
model = sprintf('%s(%i,%i,%i)%i', model, o.P(i), o.D(i), o.Q(i), o.S(i));
end
% Set random name for temporary *.spc file.
basefilename = randomstring(10);
% Write *.spc file for estimation.
fid = fopen([basefilename '.spc'], 'w');
fprintf(fid, '# File created on %s by Dynare.\n\n', datetime());
fprintf(fid, 'series {\n');
fprintf(fid, ' title = "%s"\n', o.y.name{1});
fprintf(fid, ' start = %i.%i\n', o.y.init.year, o.y.init.subperiod);
fprintf(fid, ' period = %i\n', o.y.init.freq);
fprintf(fid, ' data = %s', sprintf(data2txt(o.y.data)));
fprintf(fid, '}\n\n');
fprintf(fid, 'arima{ Model = %s }\n\n', model);
fprintf(fid, 'estimate{\n');
fprintf(fid, ' tol = %d\n', options.tol);
fprintf(fid, ' maxiter = (%i)\n', options.maxiter);
fprintf(fid, ' exact = %s\n', options.exact);
if options.outofsample
fprintf(fid, ' outofsample = yes\n');
else
fprintf(fid, ' outofsample = no\n');
end
fprintf(fid, ' print = (%s)\n', options.print);
fprintf(fid, ' save = (%s)\n', options.save);
fprintf(fid, '}\n');
fclose(fid);
% Run estimation of the ARIMA model
[status, result] = system(sprintf('%s %s', select_x13_binary(), basefilename));
% Get the content of the generated *.est file
fid = fopen(sprintf('%s.est', basefilename), 'r');
if fid<=0, error('arima:estimate: Unable to find %.est', basefilename); end
txt = textscan(fid, '%s', 'delimiter', '\n');
txt = txt{1};
fclose(fid);
% Locate the lines where the estimates are reported.
l0 = find(strncmp(txt, '$arima$estimates$:',7))+3;
l1 = find(strncmp(txt, '$variance$:',7))-1;
% Read and store the estimation results.
for l = l0:l1
linea = textscan(txt{l}, '%s %s %d %d %f %f');
switch linea{1}{1}
case 'AR'
switch linea{2}{1}
case 'Nonseasonal'
o.ar(linea{4}) = linea{5};
o.ar_std(linea{4}) = linea{6};
case 'Seasonal'
i = find(linea{3}==o.S);
o.AR{i}(linea{4}) = linea{5};
o.AR_std{i}(linea{4}) = linea{6};
otherwise
% For seasonal factors for wich the period is not matching the data frequency.
linea = textscan(txt{l}, '%s %s %d %d %d %f %f');
i = find(linea{4}==o.S);
o.AR{i}(linea{5}) = linea{6};
o.AR_std{i}(linea{5}) = linea{7};
end
case 'MA'
switch linea{2}{1}
case 'Nonseasonal'
o.ma(linea{4}) = linea{5};
o.ma_std(linea{4}) = linea{6};
case 'Seasonal'
i = find(linea{4}==o.S);
o.MA{i}(linea{4}) = linea{5};
o.MA_std{i}(linea{4}) = linea{6};
otherwise
% For seasonal factors for wich the period is not matching the data frequency.
linea = textscan(txt{l}, '%s %s %d %d %d %f %f');
i = find(linea{4}==o.S);
o.MA{i}(linea{5}) = linea{6};
o.MA_std{i}(linea{5}) = linea{7};
end
otherwise
error('arima:estimate: Unable to interpret line %i line in %s.est!', l, basefilename)
end
end
l2 = find(strncmp(txt, '$variance$:',7))+2;
% Get the value of the likelihood for the estimated parameters.
linea = textscan(txt{l2}, '%s %f');
if isequal(linea{1}{1}, 'mle')
o.estimation.mle = linea{2};
else
error('arima:estimate: Unable to find the value of the likelihood in %s.est', basefilename)
end
% Get the square root variance of the estimated innovations.
linea = textscan(txt{l2+1}, '%s %f');
if isequal(linea{1}{1}, 'se')
o.sigma = linea{2}
else
error('arima:estimate: Unable to find the estimated standard deviation of the innovations in %s.est', basefilename)
end
% Get the content of the generated *.rsd file
fid = fopen(sprintf('%s.rsd', basefilename), 'r');
if fid<=0, error('arima:estimate: Unable to find %.rsd', basefilename); end
txt = textscan(fid, '%s', 'delimiter', '\n');
txt = txt{1}(3:end);
fclose(fid);
% Read and store the estimated innovations.
data = cellfun(@(x) x{2}, cellfun(@(x) [textscan(x, '%d %f')], txt, 'UniformOutput', false), 'UniformOutput', false);
data = transpose([data{:}]);
o.estimation.e = dseries(data, o.y.init, {'innovations'}, {'\hat{\varepsilon}'});
% Get the content of the generated *.lks file
fid = fopen(sprintf('%s.lks', basefilename), 'r');
if fid<=0, error('arima:estimate: Unable to find %.lks', basefilename); end
txt = textscan(fid, '%s', 'delimiter', '\n');
fclose(fid);
linea = textscan(txt{1}{2}, '%s %d');
if isequal(linea{1}{1}, 'nobs')
o.estimation.number_of_observations = linea{2};
end
linea = textscan(txt{1}{3}, '%s %d');
if isequal(linea{1}{1}, 'nefobs')
o.estimation.effective_number_of_observations = linea{2};
end
linea = textscan(txt{1}{5}, '%s %d');
if isequal(linea{1}{1}, 'np')
o.estimation.number_of_estimated_parameters = linea{2};
end
% TODO Check the meaning of lnlkhd (versus mle)
% linea = textscan(txt{1}{8}, '%s %d');
% if isequal(linea{1}{1}, 'lnlkhd')
% o.estimation.log_likelihood = linea{2};
% end
linea = textscan(txt{1}{9}, '%s %d');
if isequal(linea{1}{1}, 'aic')
o.estimation.information_criteria.Akaike = linea{2};
end
linea = textscan(txt{1}{10}, '%s %d');
if isequal(linea{1}{1}, 'Aicc')
% Akaike with correction for finite sample size.
o.estimation.information_criteria.CorrectedAkaike = linea{2};
end
linea = textscan(txt{1}{11}, '%s %d');
if isequal(linea{1}{1}, 'hnquin')
o.estimation.information_criteria.HannanQuinn = linea{2};
end
linea = textscan(txt{1}{12}, '%s %d');
if isequal(linea{1}{1}, 'bic')
o.estimation.information_criteria.Schwarz = linea{2};
end
......@@ -19,9 +19,11 @@ dseries_src_root = strrep(which('initialize_dseries_toolbox'),'initialize_dserie
addpath([dseries_src_root '/read'])
addpath([dseries_src_root '/utilities/is'])
addpath([dseries_src_root '/utilities/str'])
addpath([dseries_src_root '/utilities/x13'])
addpath([dseries_src_root '/utilities/insert'])
addpath([dseries_src_root '/utilities/file'])
addpath([dseries_src_root '/utilities/from'])
addpath([dseries_src_root '/utilities/print'])
addpath([dseries_src_root '/utilities/variables'])
addpath([dseries_src_root '/utilities/cumulate'])
......
function b = is64bit()
% Returns true iff Matlab 64bit version is used.
% Copyright (C) 2017 Dynare Team
%
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare dseries submodule is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
[~, maxsize] = computer();
b = (maxsize > 2^31);
\ No newline at end of file
function str= data2txt(y)
% Prints a one dimensional array.
% Copyright (C) 2017 Dynare Team
%
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare dates submodule is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
str = mat2str(y);
str = strrep(str, '[', '(\n');
str = strrep(str, ';', '\n');
str = strrep(str, ']', '\n)');
\ No newline at end of file
function x13_binary = select_x13_binary()
% Copyright (C) 2017 Dynare Team
%
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare dates submodule is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
dseries_src_root = strrep(which('initialize_dseries_toolbox'),'initialize_dseries_toolbox.m','');
dseries_x13_root = sprintf('%s%s%s%s%s%s%s', dseries_src_root, '..', filesep(), 'externals', filesep(), 'x13', filesep());
if isunix()
x13_binary = sprintf('%s%s%s', dseries_x13_root, 'linux', filesep());
if is64bit()
x13_binary = sprintf('%s%s%s%s', x13_binary, '64', filesep(), 'x13');
else
x13_binary = sprintf('%s%s%s%s', x13_binary, '32', filesep(), 'x13');
end
elseif ispc()
x13_binary = sprintf('%s%s%s', dseries_x13_root, 'windows', filesep());
if is64bit()
x13_binary = sprintf('%s%s%s%s', x13_binary, '64', filesep(), 'x13.exe');
else
x13_binary = sprintf('%s%s%s%s', x13_binary, '32', filesep(), 'x13.exe');
end
else
error('X13 binary is not yet available for this plateform')
end
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment