matlab 的arima函数

classdef (Sealed) arima < internal.econ.LagIndexableTimeSeries %ARIMA Create an ARIMA model % % Syntax: % % OBJ = arima(p,D,q) % OBJ = arima(param1,val1,param2,val2,...) % % Description: % % Create an ARIMA(p,D,q) model by specifying either the degrees p, D, and % q (short-hand syntax for non-seasonal models) or a list of parameter % name-value pairs (long-hand syntax). For response process y(t) and model % innovations e(t), either syntax allows the creation of an ARIMA model of % the general form % % y(t) = c + a1*y(t-1) + ... + aP*y(t-P) + e(t) + b1*e(t-1) + ... + bQ*e(t-Q) % % with P auto-regressive terms and Q moving average terms (see notes below % for a discussion of the relationship between inputs and model degrees). % % Input Arguments (Short-Hand Syntax for Non-Seasonal Models Only): % % p - Nonnegative integer indicating the degree of the non-seasonal % auto-regressive polynomial. % % D - Nonnegative integer indicating the degree of the non-seasonal % differencing polynomial (the degree of non-seasonal integration). % % q - Nonnegative integer indicating the degree of the non-seasonal % moving average polynomial. % % Input Arguments (Parameter Name/Value Pairs): % % 'Constant' Scalar constant c of the model. If unspecified, the constant % is set to NaN. % % 'AR' A cell vector of non-seasonal auto-regressive coefficients. % When specified without corresponding lags, AR is a cell % vector of coefficients at lags 1, 2, ... to the degree of % the non-seasonal auto-regressive polynomial. When specified % along with ARLags (see below), AR is a commensurate length % cell vector of coefficients associated with the lags in % ARLags. If unspecified, AR is a cell vector of NaNs the % same length as ARLags (see below). % % 'MA' A cell vector of non-seasonal moving average coefficients. % When specified without corresponding lags, MA is a cell % vector of coefficients at lags 1, 2, ... to the degree of % the non-seasonal moving average polynomial. When specified % along with MALags (see below), MA is a commensurate length % cell vector of coefficients associated with the lags in % MALags. If unspecified, MA is a cell vector of NaNs the % same length as MALags (see below). % % 'ARLags' A vector of positive integer lags associated with the AR % coefficients. If unspecified, ARLags is a vector of integers % 1, 2, ... to the degree of the non-seasonal auto-regressive % polynomial (see AR above). % % 'MALags' A vector of positive integer lags associated with the MA % coefficients. If unspecified, MALags is a vector of integers % 1, 2, ... to the degree of the non-seasonal moving average % polynomial (see MA above). % % 'SAR' A cell vector of seasonal auto-regressive coefficients. When % specified without corresponding lags, SAR is a cell vector % of coefficients at lags 1, 2, ... to the degree of the % seasonal auto-regressive polynomial. When specified along % with SARLags (see below), SAR is a commensurate length cell % vector of coefficients associated with the lags in SARLags. % If unspecified, SAR is a cell vector of NaNs the same length % as SARLags (see below). % % 'SMA' A cell vector of seasonal moving average coefficients. When % specified without corresponding lags, SMA is a cell vector % of coefficients at lags 1, 2, ... to the degree of the % seasonal moving average polynomial. When specified along % with SMALags (see below), SMA is a commensurate length cell % vector of coefficients associated with the lags in SMALags. % If unspecified, SMA is a cell vector of NaNs the same length % as SMALags (see below). % % 'SARLags' A vector of positive integer lags associated with the SAR % coefficients. If unspecified, SARLags is a vector of integers % 1, 2, ... to the degree of the seasonal auto-regressive % polynomial (see SAR above) % % 'SMALags' A vector of positive integer lags associated with the SMA % coefficients. If unspecified, SMALags is a vector of integers % 1, 2, ... to the degree of the seasonal moving average % polynomial (see SMA above) % % 'D' A nonnegative integer indicating the degree of the non- % seasonal differencing polynomial (the degree of non-seasonal % integration). If unspecified, the default is zero. % %'Seasonality' A nonnegative integer indicating the degree of the seasonal % differencing polynomial (the degree of seasonal integration). % If unspecified, the default is zero. % % 'Beta' A vector of regression coefficients associated with a % regression component in the conditional mean. The presence % of a non-empty Beta vector allows for an ARIMAX model in % which exogenous data is included in the ARIMA model equation % shown above. % % 'Variance' A positive scalar variance of the model innovations, or a % supported conditional variance model object (e.g., a GARCH % model). If unspecified, the default is NaN. % %'Distribution' The conditional probability distribution of the innovations % process. Distribution is a string specified as 'Gaussian' % or 't', or a structure with field 'Name' which stores % the distribution 'Gaussian' or 't'. If the distribution is % 't', then the structure must also have the field 'DoF' to % store the degrees-of-freedom. % % Output Argument: % % OBJ - An ARIMA model with the following properties: % % o Distribution % o P % o D % o Q % o Constant % o AR % o MA % o SAR % o SMA % o Beta % o Variance % o Seasonality % % Notes: % % o The properties P and Q of ARIMA models do not necessarily conform to % standard Box and Jenkins notation. % % The property P is the degree of the compound autoregressive polynomial, % or the total number of lagged observations of the underlying process % necessary to initialize the auto-regressive component of the model. % P includes the effects of non-seasonal and seasonal integration % captured by the properties D and Seasonality, respectively, and the % non-seasonal and seasonal auto-regressive polynomials AR and SAR, % respectively. % % The property Q is the degree of the compound moving average polynomial, % or the total number of lagged innovations of the underlying process % necessary to initialize the moving average component of the model. Q % includes the effects of non-seasonal and seasonal moving average % polynomials MA and SMA, respectively. % % If the model has no integration and no seasonal components, only then % will the properties P and Q conform to standard Box and Jenkins % notation for an ARIMA(P,0,Q) = ARMA(P,Q) model. % % o The lags associated with the seasonal polynomials SAR and SMA are % specified in the periodicity of the observed data (as are AR and MA), % and not as multiples of the Seasonality parameter. This convention % does not conform to standard Box and Jenkins notation, yet is a more % flexible approach to incorporate multiplicative seasonal models. % % o The coefficients AR, SAR, MA, and SMA are each associated with an % underlying lag operator polynomial and subject to a near-zero % tolerance exclusion test. That is, each coefficient is compared to % the default zero tolerance 1e-12, and is included in the model only % if the magnitude is greater than 1e-12; if the coefficient magnitude % is less than or equal to 1e-12, then it is sufficiently close to zero % and excluded from the model. See LagOp for additional details. % % See also ESTIMATE, FORECAST, INFER, SIMULATE, FILTER, IMPULSE. % Copyright 1999-2011 The MathWorks, Inc. % $Revision: 1.1.6.10 $ $Date: 2012/11/15 13:42:51 $ properties (GetAccess = public, SetAccess = private, Dependent) P % Degree of the composite autoregressive polynomial Q % Degree of the composite moving average polynomial end properties (Dependent) AR % Non-seasonal autoregressive coefficients SAR % Seasonal autoregressive coefficients MA % Non-seasonal moving average coefficients SMA % Seasonal moving average coefficients end properties (Access = protected) LHS = {'AR' 'SAR'} % Published list of LHS coefficients to be reflected. RHS = {'MA' 'SMA'} % Published list of RHS coefficients NOT to be reflected. end properties (Access = public) D = 0 % Degree of non-seasonal integration Seasonality = 0 % Degree of seasonal integration Constant = NaN % Model constant Variance = NaN % Model variance Beta = [] % Regression coefficients end properties (Access = private) LagOpLHS = [] LagOpRHS = [] end methods % GET/SET methods % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function coefficients = get.AR(OBJ) coefficients = toCellArray(reflect(OBJ.LagOpLHS{1})); coefficients = coefficients(2:end); % Coefficients at lags 1, 2, ... end % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function coefficients = get.SAR(OBJ) coefficients = toCellArray(reflect(OBJ.LagOpLHS{2})); coefficients = coefficients(2:end); % Coefficients at lags 1, 2, ... end % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function coefficients = get.MA(OBJ) coefficients = toCellArray(OBJ.LagOpRHS{1}); coefficients = coefficients(2:end); % Coefficients at lags 1, 2, ... end % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function coefficients = get.SMA(OBJ) coefficients = toCellArray(OBJ.LagOpRHS{2}); coefficients = coefficients(2:end); % Coefficients at lags 1, 2, ... end % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function P = get.P(OBJ) P = OBJ.LagOpLHS{1}.Degree + OBJ.LagOpLHS{2}.Degree + OBJ.D + OBJ.Seasonality; end % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function Q = get.Q(OBJ) Q = OBJ.LagOpRHS{1}.Degree + OBJ.LagOpRHS{2}.Degree; end % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function OBJ = set.Constant(OBJ, value) validateattributes(value, {'double'}, {'scalar'}, '', 'Constant'); OBJ.Constant = value; end % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function polynomial = getLagOp(OBJ, name) %GETLAGOP Get lag operator polynomials of ARIMA models. % % Syntax: % % polynomial = getLagOp(OBJ,name) % % Description: % % Given the name of a component polynomial of the input model, return the % underlying lag operator polynomial (LagOp object). % % Input Arguments: % % OBJ - The model whose component polynomial is requested. % % name - The name of the polynomial (case-sensitive character string). % The available options are 'AR' (non-seasonal autoregressive), 'SAR' % (seasonal autoregressive), 'MA' (non-seasonal moving average), 'SMA' % (seasonal moving average), 'Integrated Non-Seasonal' (polynomial % associated with D, the degree of non-seasonal differencing), % 'Integrated Seasonal' (the polynomial associated with Seasonality, % the degree of non-seasonal differencing), 'Compound AR' (the compound % autoregressive polynomial, including all seasonal and non-seasonal % components), and 'Compound MA' (the compound moving average polynomial, % including all seasonal and non-seasonal components). % % Output Arguments: % % polynomial - A lag operator polynomial associated with the input name. switch name % Name of the component polynomial case OBJ.LHS{1} % AR (Non-Seasonal Autoregressive) polynomial = OBJ.LagOpLHS{1}; case OBJ.LHS{2} % SAR (Seasonal Autoregressive) polynomial = OBJ.LagOpLHS{2}; case OBJ.RHS{1} % MA (Non-Seasonal Moving Average) polynomial = OBJ.LagOpRHS{1}; case OBJ.RHS{2} % SMA (Seasonal Moving Average) polynomial = OBJ.LagOpRHS{2}; case 'Integrated Non-Seasonal' D = OBJ.D; if D > 0 polynomial = LagOp([1 -1]); for i = 2:D polynomial = polynomial * LagOp([1 -1]); end else polynomial = LagOp(1); end case 'Integrated Seasonal' S = OBJ.Seasonality; if S > 0 polynomial = LagOp([1 -1], 'Lags', [0 S]); else polynomial = LagOp(1); end case 'Compound AR' polynomial = OBJ.LagOpLHS{1} * OBJ.getLagOp('Integrated Seasonal') * ... OBJ.LagOpLHS{2} * OBJ.getLagOp('Integrated Non-Seasonal'); case 'Compound MA' polynomial = OBJ.LagOpRHS{1} * OBJ.LagOpRHS{2}; otherwise error(message('econ:arima:getLagOp:InvalidPolynomialReference')) end end % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function OBJ = setLagOp(OBJ, name, polynomial) %SETLAGOP Set lag operator polynomials of ARIMA models. % % Syntax: % % OBJ = setLagOp(OBJ,name,polynomial) % % Description: % % Given the name of a component polynomial of the input model, update the % underlying lag operator polynomial (LagOp object). % % Input Arguments: % % OBJ - The model whose component polynomial is updated. % % name - The name of the polynomial (case-sensitive character string). % The available options are: 'AR' (non-seasonal autoregressive), 'SAR' % (seasonal autoregressive), 'MA' (non-seasonal moving average), and % 'SMA' (seasonal moving average). % % polynomial - A lag operator polynomial associated with the input name. % % Output Arguments: % % OBJ - The model whose component polynomial is updated. switch name % Name of the component polynomial case OBJ.LHS{1} % AR polynomial OBJ.LagOpLHS{1} = polynomial; case OBJ.LHS{2} % SAR polynomial OBJ.LagOpLHS{2} = polynomial; case OBJ.RHS{1} % MA polynomial OBJ.LagOpRHS{1} = polynomial; case OBJ.RHS{2} % SMA polynomial OBJ.LagOpRHS{2} = polynomial; otherwise error(message('econ:arima:setLagOp:InvalidPolynomialReference')) end OBJ = validateModel(OBJ); end end % METHODS Block methods (Access = public) % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function OBJ = arima(varargin) %ARIMA Construct an ARIMA(P,D,Q) model if nargin == 0 % % MATLAB classes should construct a scalar object in a default state in % the absence of input arguments. % % In this case, the default constructor syntax creates an ARIMA(0,0,0) model % with a Gaussian conditional probability distribution, an undefined additive % constant, and an undefined constant-variance model. % % Therefore, the default ARIMA model represents is time series of mean-zero % unit-variance residuals. % OBJ.LagOpLHS = {LagOp(1) LagOp(1)}; OBJ.LagOpRHS = {LagOp(1) LagOp(1)}; OBJ.Constant = NaN; OBJ.Variance = NaN; return end % % Validate input parameters. % if isnumeric(varargin{1}) % Is it the short-hand syntax? % % Validate short-hand syntax, OBJ = ARIMA(P,D,Q). % parser = inputParser; parser.CaseSensitive = true; parser.addRequired('P', @(x) validateattributes(x, {'double'}, {'scalar' 'nonnegative' 'integer'}, '', 'P')); parser.addRequired('D', @(x) validateattributes(x, {'double'}, {'scalar' 'nonnegative' 'integer'}, '', 'D')); parser.addRequired('Q', @(x) validateattributes(x, {'double'}, {'scalar' 'nonnegative' 'integer'}, '', 'Q')); parser.parse(varargin{:}); P = parser.Results.P; Q = parser.Results.Q; OBJ.D = parser.Results.D; OBJ.LagOpLHS = {LagOp([1 nan(1,P)]) LagOp(1)}; OBJ.LagOpRHS = {LagOp([1 nan(1,Q)]) LagOp(1)}; elseif ischar(varargin{1}) % % Validate long-hand syntax in which the user must specify a list of % parameter name-value pairs. % OBJ = validateModel(OBJ, varargin{:}); else error(message('econ:arima:arima:InvalidInputSyntax')) end end % Constructor % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function print(OBJ, covariance) %PRINT Print ARIMA model estimation results % % Syntax: % % print(OBJ, VarCov) % % Description: % % Given an ARIMA model containing parameter estimates and corresponding % estimation error variance-covariance matrix, print a table of parameter % estimates, standard errors, and t-statistics. % % Input Arguments: % % OBJ - ARIMA model specification. % % VarCov - Estimation error variance-covariance matrix. VarCov is a square % matrix with a row and column associated with each parameter known to % the optimizer. Known parameters include all parameters estimated as % well as all coefficients held fixed throughout the optimization. % Coefficients of lag operator polynomials at lags excluded from the % model (held fixed at zero) are omitted from VarCov. % % Note: % % o The parameters known to FMINCON and included in VarCov are ordered as % follows: % % - Constant % - Non-zero AR coefficients at positive lags % - Non-zero SAR coefficients at positive lags % - Non-zero MA coefficients at positive lags % - Non-zero SMA coefficients at positive lags % - Regression coefficients (models with regression components only) % - Variance parameters (scalar for constant-variance models, vector % of additional parameters otherwise) % - Degrees-of-freedom (t distributions only) % if size(covariance,1) ~= size(covariance,2) error(message('econ:arima:print:NonSquareCovarianceMatrix')) end % % Determine which parameters are estimated and which are held fixed. % % Parameters held fixed throughout estimation are associated with all-zero % rows and columns in the error covariance matrix. % solve = sum(covariance,2) ~= 0; % TRUE = estimated, FALSE = fixed % % Get model parameters and ensure coefficients are specified. % constant = OBJ.Constant; % Model constant beta = OBJ.Beta; % Regression coefficients beta = beta(:)'; % Guarantee a column vector AR = reflect(getLagOp(OBJ, 'AR')); SAR = reflect(getLagOp(OBJ, 'SAR')); MA = getLagOp(OBJ, 'MA'); SMA = getLagOp(OBJ, 'SMA'); LagsAR = AR.Lags; % Non-zero lags of non-seasonal AR coefficients LagsSAR = SAR.Lags; % Non-zero lags of seasonal AR coefficients LagsMA = MA.Lags; % Non-zero lags of non-seasonal MA coefficients LagsSMA = SMA.Lags; % Non-zero lags of seasonal MA coefficients LagsAR = LagsAR(LagsAR > 0); % Retain only positive lags LagsSAR = LagsSAR(LagsSAR > 0); LagsMA = LagsMA(LagsMA > 0); LagsSMA = LagsSMA(LagsSMA > 0); if isempty(LagsAR) AR = []; else AR = AR.Coefficients; AR = [AR{LagsAR}]; % Non-zero AR coefficients (vector) end if isempty(LagsSAR) SAR = []; else SAR = SAR.Coefficients; SAR = [SAR{LagsSAR}]; % Non-zero SAR coefficients (vector) end if isempty(LagsMA) MA = []; else MA = MA.Coefficients; MA = [MA{LagsMA}]; % Non-zero MA coefficients (vector) end if isempty(LagsSMA) SMA = []; else SMA = SMA.Coefficients; SMA = [SMA{LagsSMA}]; % Non-zero SMA coefficients (vector) end isDistributionT = strcmpi(OBJ.Distribution.Name, 'T'); isVarianceConstant = ~isa(OBJ.Variance, 'internal.econ.LagIndexableTimeSeries'); % Is it a constant variance model? % % Determine the number of coefficients associated with the ARIMA model (i.e., % excluding any parameters found in OBJ associated with the variance model % and the degrees-of-freedom of t distributions). % nARIMA = 1 + numel(LagsAR) + numel(LagsSAR) + ... numel(LagsMA) + numel(LagsSMA) + numel(beta); % % Pack the vector of parameters associated with the ARIMA model; each element % of the parameters vector will be printed by THIS method. % % In the event a t distribution is encountered, notice that the last element % of the parameters vector is always the degrees-of-freedom regardless of % whether the variance model is constant or conditional. % if isVarianceConstant % % For constant-variance models, the number of elements in PARAMETERS will % equal nARIMA + 1 for Gaussian distributions and nARIMA + 2 for t distributions. % if isDistributionT parameters = zeros(nARIMA + 2,1); parameters(end - 1) = OBJ.Variance; parameters(end) = OBJ.Distribution.DoF; else parameters = zeros(nARIMA + 1,1); parameters(end) = OBJ.Variance; end parameters(1:nARIMA) = [constant AR SAR MA SMA beta]; if numel(parameters) ~= numel(solve) error(message('econ:arima:print:ModelCovarianceInconsistency')) end else % Conditional variance model % % For conditional-variance models, the number of elements in PARAMETERS % will be fewer than the number of elements in the SOLVE indicator because % the coefficients associated with the variance model are printed by the % print method of the variance model, and are therefore excluded from % PARAMETERS. % if isDistributionT parameters(nARIMA + 1) = OBJ.Distribution.DoF; parameters(1:nARIMA) = [constant AR SAR MA SMA beta]; else parameters = [constant AR SAR MA SMA beta]; end end % % Display summary information. % disp(' ') if isempty(OBJ.Beta) name = 'ARIMA'; else name = 'ARIMAX'; end if OBJ.Seasonality > 0 s = [' ' name '(%d,%d,%d) Model Seasonally Integrated']; else s = [' ' name '(%d,%d,%d) Model']; end if (OBJ.LagOpLHS{2}.Degree > 0) && (OBJ.LagOpRHS{2}.Degree > 0) s = [s ' with Seasonal AR(%d) and MA(%d):\n']; fprintf(s, OBJ.LagOpLHS{1}.Degree, OBJ.D, OBJ.LagOpRHS{1}.Degree, ... OBJ.LagOpLHS{2}.Degree, OBJ.LagOpRHS{2}.Degree) fprintf(' %s\n', repmat('-', 1, numel(s) - 8)) elseif OBJ.LagOpLHS{2}.Degree > 0 s = [s ' with Seasonal AR(%d):\n']; fprintf(s, OBJ.LagOpLHS{1}.Degree, OBJ.D, OBJ.LagOpRHS{1}.Degree, ... OBJ.LagOpLHS{2}.Degree) fprintf(' %s\n', repmat('-', 1, numel(s) - 8)) elseif OBJ.LagOpRHS{2}.Degree > 0 s = [s ' with Seasonal MA(%d):\n']; fprintf(s, OBJ.LagOpLHS{1}.Degree, OBJ.D, OBJ.LagOpRHS{1}.Degree, ... OBJ.LagOpRHS{2}.Degree) fprintf(' %s\n', repmat('-', 1, numel(s) - 8)) else s = [s ':\n']; fprintf(s, OBJ.LagOpLHS{1}.Degree, OBJ.D, OBJ.LagOpRHS{1}.Degree) fprintf(' %s\n', repmat('-', 1, numel(s) - 8)) end fprintf(' Conditional Probability Distribution: %s\n\n' , OBJ.Distribution.Name); header = [' Standard t ' ; ' Parameter Value Error Statistic ' ; ' ----------- ----------- ------------ -----------']; disp(header) % % Format annotation strings and display the information. Standard errors % and t-statistics held fixed by equality constraints during estimation % have zero estimation error, and are designated 'Fixed'. % Fix = ~solve; errors = sqrt(diag(covariance)); if Fix(1) fprintf(' Constant %12.6g %12s %12s\n', parameters(1), 'Fixed', 'Fixed') else tStatistic = parameters(1) / errors(1); % Compute t-statistic for this parameter fprintf(' Constant %12.6g %12.6g %12.6g\n', parameters(1), errors(1), tStatistic) end row = 2; for i = LagsAR s = [' ', repmat(' ', 1, (0 <= i) && (i < 10))]; if Fix(row) fprintf('%sAR{%d} %12.6g %12s %12s \n', s, i, parameters(row), 'Fixed', 'Fixed') else tStatistic = parameters(row) / errors(row); fprintf('%sAR{%d} %12.6g %12.6g %12.6g\n', s, i, parameters(row), errors(row), tStatistic) end row = row + 1; end for i = LagsSAR s = [' ', repmat(' ', 1, (0 <= i) && (i < 10))]; if Fix(row) fprintf('%sSAR{%d} %12.6g %12s %12s\n', s, i, parameters(row), 'Fixed', 'Fixed') else tStatistic = parameters(row) / errors(row); fprintf('%sSAR{%d} %12.6g %12.6g %12.6g\n', s, i, parameters(row), errors(row), tStatistic) end row = row + 1; end for i = LagsMA s = [' ', repmat(' ', 1, (0 <= i) && (i < 10))]; if Fix(row) fprintf('%sMA{%d} %12.6g %12s %12s \n', s, i, parameters(row), 'Fixed', 'Fixed') else tStatistic = parameters(row) / errors(row); fprintf('%sMA{%d} %12.6g %12.6g %12.6g\n', s, i, parameters(row), errors(row), tStatistic) end row = row + 1; end for i = LagsSMA s = [' ', repmat(' ', 1, (0 <= i) && (i < 10))]; if Fix(row) fprintf('%sSMA{%d} %12.6g %12s %12s\n', s, i, parameters(row), 'Fixed', 'Fixed') else tStatistic = parameters(row) / errors(row); fprintf('%sSMA{%d} %12.6g %12.6g %12.6g\n', s, i, parameters(row), errors(row), tStatistic) end row = row + 1; end for i = 1:numel(beta) s = [' ', repmat(' ', 1, (0 <= i) && (i < 10))]; if Fix(row) fprintf('%sBeta%d %12.6g %12s %12s\n', s, i, parameters(row), 'Fixed', 'Fixed') else tStatistic = parameters(row) / errors(row); fprintf('%sBeta%d %12.6g %12.6g %12.6g\n', s, i, parameters(row), errors(row), tStatistic) end row = row + 1; end if isVarianceConstant if Fix(row) fprintf(' Variance %12.6g %12s %12s\n', parameters(row), 'Fixed', 'Fixed') else tStatistic = parameters(row) / errors(row); fprintf(' Variance %12.6g %12.6g %12.6g\n', parameters(row), errors(row), tStatistic) end if isDistributionT if Fix(end) fprintf(' DoF %12.6g %12s %12s\n', parameters(end), 'Fixed', 'Fixed') else tStatistic = parameters(end) / errors(end); fprintf(' DoF %12.6g %12.6g %12.6g\n', parameters(end), errors(end), tStatistic) end end else % Conditional variance model if isDistributionT if Fix(end) fprintf(' DoF %12.6g %12s %12s\n', parameters(end), 'Fixed', 'Fixed') else tStatistic = parameters(end) / errors(end); fprintf(' DoF %12.6g %12.6g %12.6g\n', parameters(end), errors(end), tStatistic) end end disp(' ') try OBJ.Variance.print(covariance((nARIMA + 1):end,(nARIMA + 1):end)) catch exception exception.throwAsCaller(); end end end % Print Method % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function disp(OBJ) %DISP Display ARIMA models spaces = ' '; if isempty(OBJ.Beta) name = 'ARIMA'; else name = 'ARIMAX'; end if OBJ.Seasonality > 0 s = [' ' name '(%d,%d,%d) Model Seasonally Integrated']; else s = [' ' name '(%d,%d,%d) Model']; end if (OBJ.LagOpLHS{2}.Degree > 0) && (OBJ.LagOpRHS{2}.Degree > 0) s = [s ' with Seasonal AR(%d) and MA(%d):\n']; fprintf(s, OBJ.LagOpLHS{1}.Degree, OBJ.D, OBJ.LagOpRHS{1}.Degree, ... OBJ.LagOpLHS{2}.Degree, OBJ.LagOpRHS{2}.Degree) fprintf(' %s\n', repmat('-', 1, numel(s) - 8)) elseif OBJ.LagOpLHS{2}.Degree > 0 s = [s ' with Seasonal AR(%d):\n']; fprintf(s, OBJ.LagOpLHS{1}.Degree, OBJ.D, OBJ.LagOpRHS{1}.Degree, ... OBJ.LagOpLHS{2}.Degree) fprintf(' %s\n', repmat('-', 1, numel(s) - 8)) elseif OBJ.LagOpRHS{2}.Degree > 0 s = [s ' with Seasonal MA(%d):\n']; fprintf(s, OBJ.LagOpLHS{1}.Degree, OBJ.D, OBJ.LagOpRHS{1}.Degree, ... OBJ.LagOpRHS{2}.Degree) fprintf(' %s\n', repmat('-', 1, numel(s) - 8)) else s = [s ':\n']; fprintf(s, OBJ.LagOpLHS{1}.Degree, OBJ.D, OBJ.LagOpRHS{1}.Degree) fprintf(' %s\n', repmat('-', 1, numel(s) - 8)) end if strcmpi(OBJ.Distribution.Name, 'Gaussian') fprintf([spaces(1:4) 'Distribution: Name = ''%s''\n'], OBJ.Distribution.Name) else fprintf([spaces(1:4) 'Distribution: Name = ''%s'', DoF = %g\n'], OBJ.Distribution.Name, OBJ.Distribution.DoF) end fprintf([spaces(1:8) ' P: %d\n'], OBJ.P) fprintf([spaces(1:8) ' D: %d\n'], OBJ.D) fprintf([spaces(1:8) ' Q: %d\n'], OBJ.Q) fprintf([spaces(1:8) 'Constant: %g\n'], OBJ.Constant) % % Print AR coefficient information. % C = reflect(OBJ.LagOpLHS{1}); C = C.Coefficients; L = OBJ.LagOpLHS{1}.Lags; L = L(:,L > 0); N = numel(L); % Get the number of included lags. if N == 0 fprintf([spaces(1:14) 'AR: {}\n']) else format = repmat(' %g', 1, N); fprintf([spaces(1:14) 'AR: {' format(2:end) '}' ' at Lags [' format(2:end) ']\n'], C{L},L) end % % Print SAR coefficient information. % C = reflect(OBJ.LagOpLHS{2}); C = C.Coefficients; L = OBJ.LagOpLHS{2}.Lags; L = L(:,L > 0); N = numel(L); % Get the number of included lags. if N == 0 fprintf([spaces(1:13) 'SAR: {}\n']) else format = repmat(' %g', 1, N); fprintf([spaces(1:13) 'SAR: {' format(2:end) '}' ' at Lags [' format(2:end) ']\n'], C{L},L) end % % Print MA coefficient information. % C = OBJ.LagOpRHS{1}.Coefficients; L = OBJ.LagOpRHS{1}.Lags; L = L(:,L > 0); N = numel(L); % Get the number of included lags. if N == 0 fprintf([spaces(1:14) 'MA: {}\n']) else format = repmat(' %g', 1, N); fprintf([spaces(1:14) 'MA: {' format(2:end) '}' ' at Lags [' format(2:end) ']\n'], C{L},L) end % % Print SMA coefficient information. % C = OBJ.LagOpRHS{2}.Coefficients; L = OBJ.LagOpRHS{2}.Lags; L = L(:,L > 0); N = numel(L); % Get the number of included lags. if N == 0 fprintf([spaces(1:13) 'SMA: {}\n']) else format = repmat(' %g', 1, N); fprintf([spaces(1:13) 'SMA: {' format(2:end) '}' ' at Lags [' format(2:end) ']\n'], C{L},L) end % % Print information about the regression coefficients. % if ~isempty(OBJ.Beta) beta = OBJ.Beta; format = repmat(' %g', 1, numel(beta)); fprintf([spaces(1:12) 'Beta: [' format(2:end) ']\n'], beta(:).') end if OBJ.Seasonality > 0 fprintf([spaces(1:5) 'Seasonality: %d\n'], OBJ.Seasonality) end if isa(OBJ.Variance, 'double') fprintf([spaces(1:8) 'Variance: %g\n'], OBJ.Variance) elseif isa(OBJ.Variance, 'internal.econ.LagIndexableTimeSeries') fprintf([spaces(1:8) 'Variance: [%s(%d,%d) Model]\n'], upper(class(OBJ.Variance)), OBJ.Variance.P, OBJ.Variance.Q) end end % Display method end % Methods (Access = public) methods (Hidden) % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function OBJ = validateModel(OBJ, varargin) if isempty(varargin) % % Updating an existing object (no additional input arguments). % AR = OBJ.AR; MA = OBJ.MA; SAR = OBJ.SAR; SMA = OBJ.SMA; ARLags = []; % Consistent with default defined in the ELSE block below. MALags = []; SARLags = []; SMALags = []; else % % A variable-length list of parameter name-value pairs is passed when % called from the constructor. % s = struct('Name', 'Gaussian'); v = OBJ.Variance; % Ensures consistency between a contained LagIndexableTimeSeries object & its container, also a LagIndexableTimeSeries object. parser = inputParser; parser.CaseSensitive = true; parser.addParamValue('Constant' , OBJ.Constant , @(x) validateattributes(x, {'double'}, {'scalar'}, '', 'Constant')); parser.addParamValue('D' , OBJ.D , @(x) validateattributes(x, {'double'}, {'scalar' 'nonnegative' 'integer'}, '', 'D')); parser.addParamValue('AR' , {} , @(x) validateattributes(x, {'double' 'cell'}, {}, '', 'AR')); parser.addParamValue('SAR' , {} , @(x) validateattributes(x, {'double' 'cell'}, {}, '', 'SAR')); parser.addParamValue('MA' , {} , @(x) validateattributes(x, {'double' 'cell'}, {}, '', 'MA')); parser.addParamValue('SMA' , {} , @(x) validateattributes(x, {'double' 'cell'}, {}, '', 'SMA')); parser.addParamValue('Distribution', s , @(x) validateattributes(x, {'char' 'struct'}, {}, '', 'Distribution')); parser.addParamValue('ARLags' , [] , @(x) validateattributes(x, {'double'}, {'integer'}, '', 'ARLags')); parser.addParamValue('MALags' , [] , @(x) validateattributes(x, {'double'}, {'integer'}, '', 'MALags')); parser.addParamValue('SARLags' , [] , @(x) validateattributes(x, {'double'}, {'integer'}, '', 'SARLags')); parser.addParamValue('SMALags' , [] , @(x) validateattributes(x, {'double'}, {'integer'}, '', 'SMALags')); parser.addParamValue('Seasonality' , OBJ.Seasonality, @(x) validateattributes(x, {'double'}, {'scalar' 'nonnegative' 'integer'}, '', 'Seasonality')); parser.addParamValue('Variance' , v , @(x) validateattributes(x, {'double' 'garch' 'gjr' 'egarch'}, {'scalar'}, '', 'Variance')); parser.addParamValue('Beta' , OBJ.Beta , @(x) validateattributes(x, {'double'}, {}, '', 'Beta')); parser.parse(varargin{:}); AR = parser.Results.AR; SAR = parser.Results.SAR; MA = parser.Results.MA; SMA = parser.Results.SMA; ARLags = parser.Results.ARLags; MALags = parser.Results.MALags; SARLags = parser.Results.SARLags; SMALags = parser.Results.SMALags; OBJ.Constant = parser.Results.Constant; OBJ.D = parser.Results.D; OBJ.Distribution = parser.Results.Distribution; OBJ.Seasonality = parser.Results.Seasonality; OBJ.Variance = parser.Results.Variance; OBJ.Beta = parser.Results.Beta; if ~isempty(ARLags) && ~isvector(ARLags) error(message('econ:arima:validateModel:NonVectorARLags')) end if ~isempty(SARLags) && ~isvector(SARLags) error(message('econ:arima:validateModel:NonVectorSARLags')) end if ~isempty(MALags) && ~isvector(MALags) error(message('econ:arima:validateModel:NonVectorMALags')) end if ~isempty(SMALags) && ~isvector(SMALags) error(message('econ:arima:validateModel:NonVectorSMALags')) end end % % Additional checking for regression coefficients. % if ~isempty(OBJ.Beta) if ~isvector(OBJ.Beta) error(message('econ:arima:validateModel:NonVectorBeta')) else beta = OBJ.Beta; OBJ.Beta = beta(:).'; % Ensure a row vector end end % % Check the variance model. % if ~any(strcmp(class(OBJ.Variance), {'double' 'garch' 'gjr' 'egarch'})) error(message('econ:arima:validateModel:InvalidVarianceType')) end if isa(OBJ.Variance, 'double') && (OBJ.Variance <= 0) % Check constant-variance models error(message('econ:arima:validateModel:NonPositiveVariance')) end if isa(OBJ.Variance, 'internal.econ.LagIndexableTimeSeries') && (OBJ.Variance.Offset ~= 0) % Check conditional variance models error(message('econ:arima:validateModel:NonZeroVarianceOffset')) end % % Enforce consistency between probability distributions. % if isa(OBJ.Variance, 'internal.econ.LagIndexableTimeSeries') % % The ARIMA model contains a conditional variance model, so set the % distribution of the variance model to that of the ARIMA model if explicitly % specified by the user, allowing the outer, or container, ARIMA model to % take precedence. % % However, if the user did NOT specify a distribution for the ARIMA model, % then allow the distribution of the contained variance model to migrate to % the ARIMA model. % if any(strcmpi('Distribution', varargin(1:2:end))) % Did the user specify a distribution? OBJ.Variance.Distribution = OBJ.Distribution; else OBJ.Distribution = OBJ.Variance.Distribution; end end % % Check and assign model coefficients and corresponding lags. % % In the following, if the user specified coefficients without lags, then % derive lags from the coefficients. Similarly, if the user specified lags % without coefficients, then assign NaNs to coefficients consistent with % the number of lags. % % % Additional checking for AR coefficients. % if isnumeric(AR) AR = num2cell(AR); end if length(AR) == numel(AR) AR = AR(:).'; else error(message('econ:arima:validateModel:NonVectorARCoefficients')) end for i = 1:numel(AR) if ~isscalar(AR{i}) error(message('econ:arima:validateModel:NonScalarARCoefficients')) end end if ~any(strcmpi('ARLags', varargin(1:2:end))) ARLags = 0:numel(AR); else if ~any(strcmpi('AR', varargin(1:2:end))) AR = num2cell(nan(1,numel(ARLags))); end if numel(ARLags) ~= numel(AR) error(message('econ:arima:validateModel:InconsistentARSpec')) end Lags = unique(ARLags, 'first'); if any(ARLags == 0) || (numel(ARLags) ~= numel(Lags)) error(message('econ:arima:validateModel:InvalidARLags')) end ARLags = [0 ARLags]; end % % Additional checking for SAR coefficients. % if isnumeric(SAR) SAR = num2cell(SAR); end if length(SAR) == numel(SAR) SAR = SAR(:).'; else error(message('econ:arima:validateModel:NonVectorSARCoefficients')) end for i = 1:numel(SAR) if ~isscalar(SAR{i}) error(message('econ:arima:validateModel:NonScalarSARCoefficients')) end end if ~any(strcmpi('SARLags', varargin(1:2:end))) SARLags = 0:numel(SAR); else if ~any(strcmpi('SAR', varargin(1:2:end))) SAR = num2cell(nan(1,numel(SARLags))); end if numel(SARLags) ~= numel(SAR) error(message('econ:arima:validateModel:InconsistentSARSpec')) end Lags = unique(SARLags, 'first'); if any(SARLags == 0) || (numel(SARLags) ~= numel(Lags)) error(message('econ:arima:validateModel:InvalidSARLags')) end SARLags = [0 SARLags]; end % % Additional checking for MA coefficients. % if isnumeric(MA) MA = num2cell(MA); end if length(MA) == numel(MA) MA = MA(:).'; else error(message('econ:arima:validateModel:NonVectorMACoefficients')) end for i = 1:numel(MA) if ~isscalar(MA{i}) error(message('econ:arima:validateModel:NonScalarMACoefficients')) end end if ~any(strcmpi('MALags', varargin(1:2:end))) MALags = 0:numel(MA); else if ~any(strcmpi('MA', varargin(1:2:end))) MA = num2cell(nan(1,numel(MALags))); end if numel(MALags) ~= numel(MA) error(message('econ:arima:validateModel:InconsistentMASpec')) end Lags = unique(MALags); if any(MALags == 0) || (numel(MALags) ~= numel(Lags)) error(message('econ:arima:validateModel:InvalidMALags')) end MALags = [0 MALags]; end % % Additional checking for SMA coefficients. % if isnumeric(SMA) SMA = num2cell(SMA); end if length(SMA) == numel(SMA) SMA = SMA(:).'; else error(message('econ:arima:validateModel:NonVectorSMACoefficients')) end for i = 1:numel(SMA) if ~isscalar(SMA{i}) error(message('econ:arima:validateModel:NonScalarSMACoefficients')) end end if ~any(strcmpi('SMALags', varargin(1:2:end))) SMALags = 0:numel(SMA); else if ~any(strcmpi('SMA', varargin(1:2:end))) SMA = num2cell(nan(1,numel(SMALags))); end if numel(SMALags) ~= numel(SMA) error(message('econ:arima:validateModel:InconsistentSMASpec')) end Lags = unique(SMALags); if any(SMALags == 0) || (numel(SMALags) ~= numel(Lags)) error(message('econ:arima:validateModel:InvalidSMALags')) end SMALags = [0 SMALags]; end % % Now create the underlying Lag Operator Polynomials. % if isempty(OBJ.LagOpLHS) OBJ.LagOpLHS{1} = LagOp([1 cellfun(@uminus, AR, 'uniformoutput', false)], ... 'Lags', ARLags); OBJ.LagOpLHS{2} = LagOp([1 cellfun(@uminus, SAR, 'uniformoutput', false)], ... 'Lags', SARLags); end if isempty(OBJ.LagOpRHS) OBJ.LagOpRHS{1} = LagOp([1 MA], 'Lags', MALags); OBJ.LagOpRHS{2} = LagOp([1 SMA], 'Lags', SMALags); end if OBJ.LagOpLHS{1}.Coefficients{0} ~= 1 error(message('econ:arima:validateModel:InvalidAR0Coefficient')) end if OBJ.LagOpLHS{2}.Coefficients{0} ~= 1 error(message('econ:arima:validateModel:InvalidSAR0Coefficient')) end if OBJ.LagOpRHS{1}.Coefficients{0} ~= 1 error(message('econ:arima:validateModel:InvalidMA0Coefficient')) end if OBJ.LagOpRHS{2}.Coefficients{0} ~= 1 error(message('econ:arima:validateModel:InvalidSMA0Coefficient')) end % % Indicate stationarity/invertibility constraints for coefficients. % if all(~isnan(cell2mat(AR))) && ~isempty(cell2mat(AR)) && ~isStable(OBJ.LagOpLHS{1}) error(message('econ:arima:validateModel:UnstableAR')) end if all(~isnan(cell2mat(SAR))) && ~isempty(cell2mat(SAR)) && ~isStable(OBJ.LagOpLHS{2}) error(message('econ:arima:validateModel:UnstableSAR')) end if all(~isnan(cell2mat(MA))) && ~isempty(cell2mat(MA)) && ~isStable(OBJ.LagOpRHS{1}) error(message('econ:arima:validateModel:UnstableMA')) end if all(~isnan(cell2mat(SMA))) && ~isempty(cell2mat(SMA)) && ~isStable(OBJ.LagOpRHS{2}) error(message('econ:arima:validateModel:UnstableSMA')) end end % Validate Model end % Methods (Hidden) methods (Static, Hidden) % Error Checking Utilities % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function [nLogL, nLogLikelihoods, E] = nLogLikeGaussian(coefficients, Y, X, E, V, LagsAR, LagsMA, numPaths, maxPQ, T) % % Infer the residuals. % I = ones(numPaths,1); if isempty(X) for t = (maxPQ + 1):T data = [I Y(:,t - LagsAR) E(:,t - LagsMA)]; E(:,t) = data * coefficients; end else arimaCoefficients = coefficients(1:(end - size(X,1))); regressCoefficients = coefficients((end - size(X,1) + 1):end)'; for t = (maxPQ + 1):T data = [I Y(:,t - LagsAR) E(:,t - LagsMA)]; E(:,t) = data * arimaCoefficients - regressCoefficients * X(:,t); end end % % Evaluate the Gaussian negative log-likelihood objective. % nLogLikelihoods = 0.5 * (log(2*pi*V(:,(maxPQ + 1):T)) + ((E(:,(maxPQ + 1):T).^2)./V(:,(maxPQ + 1):T))); nLogL = sum(nLogLikelihoods, 2); % % Trap any degenerate negative log-likelihood values. % i = isnan(nLogL) | (imag(nLogL) ~= 0) | isinf(nLogL); if any(i) nLogL(i) = 1e20; end end % Gaussian negative log-likelihood function % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % function [nLogL, nLogLikelihoods, E] = nLogLikeT(coefficients, Y, X, E, V, LagsAR, LagsMA, numPaths, maxPQ, T) % % Infer the residuals. % I = ones(numPaths,1); if isempty(X) arimaCoefficients = coefficients(1:(end - 1)); for t = (maxPQ + 1):T data = [I Y(:,t - LagsAR) E(:,t - LagsMA)]; E(:,t) = data * arimaCoefficients; end else arimaCoefficients = coefficients(1:(end - size(X,1) - 1)); regressCoefficients = coefficients((end - size(X,1)):(end - 1))'; for t = (maxPQ + 1):T data = [I Y(:,t - LagsAR) E(:,t - LagsMA)]; E(:,t) = data * arimaCoefficients - regressCoefficients * X(:,t); end end % % Evaluate standardized Student's t negative log-likelihood function. % DoF = coefficients(end); if DoF <= 200 % % Standardized t-distributed residuals. % nLogLikelihoods = 0.5 * (log(V(:,(maxPQ + 1):T)) + (DoF + 1) * log(1 + (E(:,(maxPQ + 1):T).^2)./(V(:,(maxPQ + 1):T) * (DoF - 2)))); nLogLikelihoods = nLogLikelihoods - log(gamma((DoF + 1)/2) / (gamma(DoF/2) * sqrt(pi * (DoF - 2)))); nLogL = sum(nLogLikelihoods, 2); else % % Gaussian residuals. % nLogLikelihoods = 0.5 * (log(2*pi*V(:,(maxPQ + 1):T)) + ((E(:,(maxPQ + 1):T).^2)./V(:,(maxPQ + 1):T))); nLogL = sum(nLogLikelihoods, 2); end % % Trap any degenerate negative log-likelihood values. % i = isnan(nLogL) | (imag(nLogL) ~= 0) | isinf(nLogL); if any(i) nLogL(i) = 1e20; end end % t-distribution negative log-likelihood function end % Methods (Static, Hidden) end % Class definition