From 6997e0a4a65763ef1af6cf9bdcf4864ca35f6a7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 8 Mar 2019 14:34:41 +0100 Subject: [PATCH] Added the possibility to have exogenous variables in the optimal part of PAC. Works with iterative ols and nls. --- matlab/+pac/+estimate/init.m | 3 + matlab/+pac/+estimate/iterative_ols.m | 103 ++++++++++--- tests/Makefile.am | 1 + tests/pac/trend-component-25/clean | 8 + tests/pac/trend-component-25/example.mod | 185 +++++++++++++++++++++++ 5 files changed, 276 insertions(+), 24 deletions(-) create mode 100755 tests/pac/trend-component-25/clean create mode 100644 tests/pac/trend-component-25/example.mod diff --git a/matlab/+pac/+estimate/init.m b/matlab/+pac/+estimate/init.m index b6caf2561..9b2dbb672 100644 --- a/matlab/+pac/+estimate/init.m +++ b/matlab/+pac/+estimate/init.m @@ -71,6 +71,9 @@ stack = dbstack; ipnames__ = ipnames_; % The user provided order. if isequal(stack(2).name, 'iterative_ols') ipnames_ = [M_.pac.(pacmodl).equations.(eqtag).ec.params; M_.pac.(pacmodl).equations.(eqtag).ar.params']; + if isfield(M_.pac.(pacmodl).equations.(eqtag), 'optim_additive') + ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).optim_additive.params(~isnan(M_.pac.(pacmodl).equations.(eqtag).optim_additive.params))']; + end if isfield(M_.pac.(pacmodl).equations.(eqtag), 'additive') ipnames_ = [ipnames_; M_.pac.(pacmodl).equations.(eqtag).additive.params(~isnan(M_.pac.(pacmodl).equations.(eqtag).additive.params))']; end diff --git a/matlab/+pac/+estimate/iterative_ols.m b/matlab/+pac/+estimate/iterative_ols.m index 092a6366c..d4bce328b 100644 --- a/matlab/+pac/+estimate/iterative_ols.m +++ b/matlab/+pac/+estimate/iterative_ols.m @@ -57,6 +57,13 @@ else is_exogenous_variables = false; end +% Set flag for models with exogenous variables (in the optimizing agents part) +if isfield(M_.pac.(pacmodl).equations.(eqtag), 'optim_additive') + is_optim_exogenous_variables = length(M_.pac.(pacmodl).equations.(eqtag).optim_additive.vars)>0; +else + is_optim_exogenous_variables = false; +end + if is_non_optimizing_agents non_optimizing_behaviour = M_.pac.(pacmodl).equations.(eqtag).non_optimizing_behaviour; non_optimizing_behaviour_params = NaN(length(non_optimizing_behaviour.params), 1); @@ -100,14 +107,17 @@ if is_exogenous_variables additive.scaling_factor(residual_jd) = []; additive.estimation = ismember(additive.params, ipnames_); else - additive.params = []; - additive.vars = []; - additive.isendo = []; - additive.lags = []; - additive.scaling_factor = []; - additive.estimation = []; + additive = struct('params', [], 'vars', [], 'isendo', [], 'lags', [], 'scaling_factor', [], 'estimation', []); end +if is_optim_exogenous_variables + optim_additive = M_.pac.(pacmodl).equations.(eqtag).optim_additive; + optim_additive.estimation = ismember(optim_additive.params, ipnames_); +else + optim_additive = struct('params', [], 'vars', [], 'isendo', [], 'lags', [], 'scaling_factor', [], 'estimation', []); +end + + % Build PAC expectation matrix expression. dataForPACExpectation0 = dseries(); listofvariables0 = {}; @@ -207,6 +217,43 @@ else dataForExogenousVariables_ = dseries(); end +% Build data for exogenous variables (in the optimizing behaviour term). +if is_optim_exogenous_variables + listofvariables4 = {}; j = 0; + dataForOptimExogenousVariables = dseries(); % Estimated parameters + dataForOptimExogenousVariables_ = 0; % Calibrated parameters + is_any_calibrated_parameter_optim_x = false; + is_any_estimated_parameter_optim_x = false; + for i=1:length(optim_additive.vars) + if optim_additive.isendo(i) + variable = M_.endo_names{optim_additive.vars(i)}; + else + variable = M_.exo_names{optim_additive.vars(i)}; + end + if optim_additive.estimation(i) + j = j+1; + is_any_estimated_parameter_optim_x = true; + listofvariables4{j} = variable; + dataForOptimExogenousVariables = [dataForOptimExogenousVariables, optim_additive.scaling_factor(i)*data{variable}.lag(optim_additive.lags(i))]; + else + is_any_calibrated_parameter_optim_x = true; + tmp = data{variable}.lag(optim_additive.lags(i)).data; + if ~isnan(optim_additive.params(i)) + tmp = M_.params(optim_additive.params(i))*tmp; + end + tmp = optim_additive.scaling_factor(i)*tmp; + dataForOptimExogenousVariables_ = dataForOptimExogenousVariables_+tmp; + end + end + if is_any_calibrated_parameter_optim_x + dataForOptimExogenousVariables_ = dseries(dataForOptimExogenousVariables_, data.dates(1), 'exogenous_variables_associated_with_calibrated_parameters'); + end +else + listofvariables4 = {}; + dataForOptimExogenousVariables = dseries(); + dataForOptimExogenousVariables_ = dseries(); +end + % Reorder ec.vars locally if necessary. Second variable must be the % endogenous variable, while the first must be the associated trend. if M_.pac.(pacmodl).equations.(eqtag).ec.isendo(2) @@ -232,6 +279,10 @@ end XDATA = DataForOLS{listofvariables3{:}}(range).data; +if is_optim_exogenous_variables && is_any_estimated_parameter_optim_x + XDATA = [XDATA, dataForOptimExogenousVariables{listofvariables4{:}}(range).data]; +end + if is_exogenous_variables && is_any_estimated_parameter_x XDATA = [XDATA, dataForExogenousVariables{listofvariables2{:}}(range).data]; end @@ -245,30 +296,25 @@ else end % Get indices in params0 for EC and AR parameters -[~, params_id_1] = setdiff(ipnames_, [share_of_optimizing_agents_index, additive.params]); +[~, params_id_1] = setdiff(ipnames_, [share_of_optimizing_agents_index, optim_additive.params, additive.params, ]); -% Get indices in params0 for other parameters (optimizing agents -% share plus parameters related to exogenous variables). +% Get indices in params0 for EC and AR parameters plus parameters related to exogenous variables in the optimal part. +[~, params_id_5] = setdiff(ipnames_, [share_of_optimizing_agents_index, additive.params]); + +% Get indices in params0 for other parameters (optimizing agents share plus parameters related to exogenous variables). [~, params_id_2] = setdiff(1:length(ipnames_), params_id_1); -% Get indices in params0 for the parameters associated to the -% exogenous variables. +% Get indices in params0 for the parameters associated to the exogenous variables. params_id_3 = setdiff(params_id_2, params_id_0); +[~, params_id_3_o] = ismember(optim_additive.params, ipnames_); +params_id_3_no = setdiff(params_id_3, params_id_3_o); -% Get indices in params0 for EC/AR parameters and parameters -% associated to the exogenous variables (if any). +% Get indices in params0 for EC/AR parameters and parameters associated to the exogenous variables (if any). params_id_4 = [params_id_1; params_id_3]; -[~, params_id_1] = setdiff(ipnames_, [share_of_optimizing_agents_index, additive.params]); - - -% Get values for EC-AR parameters plus the parameters associated to -% the exogenous variables (if any). +% Get values for EC-AR parameters plus the parameters associated to the exogenous variables (if any). params0_ = params0([params_id_1; params_id_3]); -% Get values for parameters associated to the exogenous variables. -params0__ = params0(params_id_3); - % Get value of the share of optimizing agents. if estimate_non_optimizing_agents_share share_of_optimizing_agents = params0(params_id_0); @@ -301,11 +347,17 @@ while noconvergence if is_exogenous_variables && is_any_calibrated_parameter_x YDATA = YDATA-dataForExogenousVariables_(range).data; end + if is_optim_exogenous_variables && is_any_calibrated_parameter_optim_x + YDATA = YDATA-share_of_optimizing_agents*dataForOptimExogenousVariables_(range).data; + end % Run OLS to estimate PAC parameters (autoregressive parameters and error correction parameter). params1_ = (XDATA\YDATA); - params1_(1:length(params_id_1)) = params1_(1:length(params_id_1))/share_of_optimizing_agents; + params1_(1:length(params_id_5)) = params1_(1:length(params_id_5))/share_of_optimizing_agents; % Compute residuals and sum of squareed residuals. - r = YDATA-XDATA*(params1_*share_of_optimizing_agents); + r = YDATA-XDATA(:,1:length(params_id_5))*params1_(1:length(params_id_5))*share_of_optimizing_agents; + if is_optim_exogenous_variables && is_any_estimated_parameter_optim_x + r = r-XDATA(:,length(params_id_5)+1:end)*params1_(length(params_id_5)+1:end); + end ssr = r'*r; % Update convergence dummy variable and display iteration info. noconvergence = max(abs(params0_-params1_))>1e-6; @@ -329,7 +381,10 @@ while noconvergence end end % Set vector for regressor. - ZDATA = XDATA(:,1:length(params_id_1))*params0_(1:length(params_id_1))+PacExpectations-dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params; + ZDATA = XDATA(:,1:length(params_id_5))*params0_(1:length(params_id_5))+PacExpectations-dataForNonOptimizingBehaviour(range).data*non_optimizing_behaviour_params; + if is_optim_exogenous_variables && is_any_calibrated_parameter_optim_x + ZDATA = ZDATA+dataForOptimExogenousVariables_(range).data; + end % Update the (estimated) share of optimizing agents by running OLS share_of_optimizing_agents = (ZDATA\YDATA); % Force the share of optimizing agents to be in [0,1]. diff --git a/tests/Makefile.am b/tests/Makefile.am index 04b85a154..e75849ae7 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -432,6 +432,7 @@ MODFILES = \ pac/trend-component-22/example.mod \ pac/trend-component-23/example.mod \ pac/trend-component-24/example.mod \ + pac/trend-component-25/example.mod \ write/example1.mod \ ecb/backward-models/irf/solow_1.mod \ ecb/backward-models/irf/solow_2.mod \ diff --git a/tests/pac/trend-component-25/clean b/tests/pac/trend-component-25/clean new file mode 100755 index 000000000..cf3492a28 --- /dev/null +++ b/tests/pac/trend-component-25/clean @@ -0,0 +1,8 @@ +#!/bin/sh + +rm -rf example +rm -rf +example +rm -f example.log +rm -f *.mat +rm -f *.m +rm -f *.dat diff --git a/tests/pac/trend-component-25/example.mod b/tests/pac/trend-component-25/example.mod new file mode 100644 index 000000000..bf609cd22 --- /dev/null +++ b/tests/pac/trend-component-25/example.mod @@ -0,0 +1,185 @@ +// --+ options: json=compute, stochastic +-- + +var x1 x2 x1bar x2bar z y x u v s; + +varexo ex1 ex2 ex1bar ex2bar ez ey ex eu ev es; + +parameters + rho_1 rho_2 rho_3 rho_4 + a_x1_0 a_x1_1 a_x1_2 a_x1_x2_1 a_x1_x2_2 + a_x2_0 a_x2_1 a_x2_2 a_x2_x1_1 a_x2_x1_2 + e_c_m c_z_1 c_z_2 c_z_dx2 c_z_u c_z_dv c_z_s beta + lambda; + +rho_1 = .9; +rho_2 = -.2; +rho_3 = .4; +rho_4 = -.3; + + +a_x1_0 = -.9; +a_x1_1 = .4; +a_x1_2 = .3; +a_x1_x2_1 = .1; +a_x1_x2_2 = .2; + +a_x2_0 = -.9; +a_x2_1 = .2; +a_x2_2 = -.1; +a_x2_x1_1 = -.1; +a_x2_x1_2 = .2; + +beta = .2; +e_c_m = .5; +c_z_1 = .2; +c_z_2 = -.1; +c_z_dx2 = .3; +c_z_u = .3; +c_z_dv = .4; +c_z_s = -.2; + +lambda = 0.5; // Share of optimizing agents. + +trend_component_model(model_name=toto, eqtags=['eq:x1', 'eq:x2', 'eq:x1bar', 'eq:x2bar'], targets=['eq:x1bar', 'eq:x2bar']); + +pac_model(auxiliary_model_name=toto, discount=beta, model_name=pacman); + +model; + +[name='eq:u'] +s = .3*s(-1) - .1*s(-2) + es; + +[name='eq:diff(v)'] +diff(v) = .5*diff(v(-1)) + ev; + +[name='eq:u'] +u = .5*u(-1) - .2*u(-2) + eu; + +[name='eq:y'] +y = rho_1*y(-1) + rho_2*y(-2) + ey; + +[name='eq:x'] +x = rho_3*x(-1) + rho_4*x(-2) + ex; + +[name='eq:x1'] +diff(x1) = a_x1_0*(x1(-1)-x1bar(-1)) + a_x1_1*diff(x1(-1)) + a_x1_2*diff(x1(-2)) + a_x1_x2_1*diff(x2(-1)) + a_x1_x2_2*diff(x2(-2)) + ex1; + +[name='eq:x2'] +diff(x2) = a_x2_0*(x2(-1)-x2bar(-1)) + a_x2_1*diff(x1(-1)) + a_x2_2*diff(x1(-2)) + a_x2_x1_1*diff(x2(-1)) + a_x2_x1_2*diff(x2(-2)) + ex2; + +[name='eq:x1bar'] +x1bar = x1bar(-1) + ex1bar; + +[name='eq:x2bar'] +x2bar = x2bar(-1) + ex2bar; + +[name='zpac'] +diff(z) = lambda*(e_c_m*(x1(-1)-z(-1)) + c_z_1*diff(z(-1)) + c_z_2*diff(z(-2)) + pac_expectation(pacman) + c_z_s*s + c_z_dv*diff(v) ) + (1-lambda)*( y + x) + c_z_u*u + c_z_dx2*diff(x2) + ez; + +end; + +shocks; + var ex1 = 1.0; + var ex2 = 1.0; + var ex1bar = 1.0; + var ex2bar = 1.0; + var ez = 1.0; + var ey = 0.1; + var ex = 0.1; + var eu = 0.05; + var ev = 0.05; + var es = 0.07; +end; + +// Initialize the PAC model (build the Companion VAR representation for the auxiliary model). +pac.initialize('pacman'); + +// Update the parameters of the PAC expectation model (h0 and h1 vectors). +pac.update.expectation('pacman'); + +// Set initial conditions to zero. Please use more sensible values if any... +initialconditions = dseries(zeros(10, M_.endo_nbr+M_.exo_nbr), 2000Q1, vertcat(M_.endo_names,M_.exo_names)); + +// Simulate the model for 500 periods +TrueData = simul_backward_model(initialconditions, 5000); + +// Define a structure describing the parameters to be estimated (with initial conditions). +clear eparams +eparams.e_c_m = .9; +eparams.c_z_1 = .5; +eparams.c_z_2 = .2; +eparams.c_z_s = .9; +eparams.c_z_dx2 = .1; +eparams.lambda = .7; + +// Define the dataset used for estimation +edata = TrueData; +edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez'); +pac.estimate.iterative_ols('zpac', eparams, edata, 2005Q1:2005Q1+4000); + +e_c_m_iterative_ols = M_.params(strmatch('e_c_m', M_.param_names, 'exact')); +c_z_1_iterative_ols = M_.params(strmatch('c_z_1', M_.param_names, 'exact')); +c_z_2_iterative_ols = M_.params(strmatch('c_z_2', M_.param_names, 'exact')); +c_z_s_iterative_ols = M_.params(strmatch('c_z_s', M_.param_names, 'exact')); +c_z_dx2_iterative_ols = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact')); +lambda_iterative_ols = M_.params(strmatch('lambda', M_.param_names, 'exact')); + +disp(sprintf('Estimate of e_c_m: %f', e_c_m_iterative_ols)) +disp(sprintf('Estimate of c_z_1: %f', c_z_1_iterative_ols)) +disp(sprintf('Estimate of c_z_2: %f', c_z_2_iterative_ols)) +disp(sprintf('Estimate of c_z_s: %f', c_z_s_iterative_ols)) +disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_iterative_ols)) +disp(sprintf('Estimate of lambda: %f', lambda_iterative_ols)) + +skipline(2) + +// Define a structure describing the parameters to be estimated (with initial conditions). +clear eparams +eparams.e_c_m = .9; +eparams.c_z_1 = .5; +eparams.c_z_2 = .2; +eparams.c_z_s = .9; +eparams.c_z_dx2 = .1; +eparams.lambda = .7; + +edata = TrueData; +edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez'); +pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon'); + +e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact')); +c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact')); +c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact')); +c_z_s_nls = M_.params(strmatch('c_z_s', M_.param_names, 'exact')); +c_z_dx2_nls = M_.params(strmatch('c_z_dx2', M_.param_names, 'exact')); +lambda_nls = M_.params(strmatch('lambda', M_.param_names, 'exact')); + +disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls)) +disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls)) +disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls)) +disp(sprintf('Estimate of c_z_s: %f', c_z_s_nls)) +disp(sprintf('Estimate of c_z_dx2: %f', c_z_dx2_nls)) +disp(sprintf('Estimate of lambda: %f', lambda_nls)) + +if abs(e_c_m_nls-e_c_m_iterative_ols)>.01 + error('Iterative OLS and NLS do not provide consistent estimates (e_c_m)') +end + +if abs(c_z_1_nls-c_z_1_iterative_ols)>.01 + error('Iterative OLS and NLS do not provide consistent estimates (c_z_1)') +end + +if abs(c_z_2_nls-c_z_2_iterative_ols)>.01 + error('Iterative OLS and NLS do not provide consistent estimates (c_z_2)') +end + +if abs(c_z_s_nls-c_z_s_iterative_ols)>.01 + error('Iterative OLS and NLS do not provide consistent estimates (c_z_s)') +end + +if abs(c_z_dx2_nls-c_z_dx2_iterative_ols)>.01 + error('Iterative OLS and NLS do not provide consistent estimates (c_z_dx2)') +end + +if abs(lambda_nls-lambda_iterative_ols)>.01 + error('Iterative OLS and NLS do not provide consistent estimates (lambda)') +end