minimizing the distance between the data and model
Show older comments
Hello, I need your help to understand the codes in below. I am trying to calibrate omega, sigma_a,sigma_g by minimizing the distance of a set of second moments with respect to the data(GMM). First, could you please tell me whether the codes in below is the right way to do that? Second, It gives me the error and do not understand the error. It says
quation solved.
fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
exitflag =
1
exitflag =
0
MOMS_hp_model (results in the file)
0.00000000000000 + 1338.61275639433i
0.00000000000000 - 0.00898886119437529i
0.00000000000000 - 0.00744518649879554i
0.00000000000000 - 0.00372482891374872i
-24.7058560639534 + 0.00000000000000i
0.00000000000000 - 0.152486798477025i
0.00000000000000 - 0.371917868654176i
0.00000000000000 - 0.0774977355076163i
0.00000000000000 - 0.00295674375382621i
0.00000000000000 - 0.312622813319814i
0.00000000000000 - 0.0759160590140773i
0.00000000000000 - 0.00518807239067139i
0.00000000000000 + 0.0881403577726612i
0.00000000000000 - 0.0126060504437991i
0.00000000000000 - 0.0399092709720994i
0.00000000000000 - 0.0541691455171444i
0.00000000000000 - 0.0734824897003621i
0.00000000000000 - 0.0910792222960299i
0.00000000000000 + 0.118458447462384i
0.00000000000000 + 0.140575614040965i
0.00000000000000 + 0.149113730723083i
0.00000000000000 + 0.148146721050559i
Here is the m file:
clear all clc
warning off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setting Empirical Moments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MOMS_hp_data =...
[2.32;...
1.31;...
3.79;...
0.78;...
0.81;...
0.92;...
0.95;...
-0.82;...
0.42;...
0.54;...
0.56;...
1.07;...
0.74;...
-0.45];
MOMS_hp_data_corrs_yF_hI = ...
[-0.2009 0.0707 0.3322;...
-0.3095 -0.0481 0.22;...
-0.4064 -0.1607 0.1067;...
-0.5245 -0.3057 -0.049;...
-0.6392 -0.4562 -0.2243;...
-0.6987 -0.5356 -0.3196;...
-0.669 -0.4927 -0.264;...
-0.66 -0.4784 -0.2441;...
-0.6368 -0.4449 -0.2011];
% Setting entire vector of moments for quadratic function
MOMS_hp_data_all = [MOMS_hp_data;...
MOMS_hp_data_corrs_yF_hI(1:4,2);...
MOMS_hp_data_corrs_yF_hI(6:9,2)];
% Grid for omega between (0,1) with equally spaced grids (except few cases
% to avoid numerical indeterminacy of the solution)
omega_grid = [0.05, 0.10, 0.15, 0.201, 0.251, 0.30, 0.349, 0.399, 0.449, 0.506, 0.551, 0.601, 0.651, 0.70, 0.752, 0.802, 0.853, 0.901, 0.95];
[x,size_omega] = size(omega_grid);
% Grid for sigma_a
sigma_a_first = 0.0001; sigma_a_last = 0.01;
sigma_a_num_grids = 1;
sigma_a_increase = (sigma_a_last-sigma_a_first)/sigma_a_num_grids;
sigma_a_grid = sigma_a_first:sigma_a_increase:sigma_a_last;
[x,size_sigma_a] = size(sigma_a_grid);
% Grid for sigma_g
sigma_g_first = 0.0001; sigma_g_last = 0.01;
sigma_g_num_grids = 1;
sigma_g_increase = (sigma_g_last-sigma_g_first)/sigma_g_num_grids;
sigma_g_grid = sigma_g_first:sigma_g_increase:sigma_g_last;
[x,size_sigma_g] = size(sigma_g_grid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calibration Benchmark Case
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calibrating some of the other key parameter values. Same calibration as
% in the benchmark case
DELTA_F = 0.05;
DELTA_I = 0.05;
TFPF_to_TFPI = 2.1901;
ALPHA_I = 0.2;
MU = 1.006;
elasticity_F_I = 8;
E = 1-(1/elasticity_F_I);
GAMMA = (MU^(-0.35+ALPHA_I)/TFPF_to_TFPI)^(1/ALPHA_I); % This calibration changes relative to the benchmark case. See appendix.
TAUN = .1142;
TAUY = .07223;
RHOA = 0.94;
RHOG = 0.72;
CHIF = 0;
CHII=0;
%%%%%%%%%%%%%%%%%%%%%%%%%
% Grid - moment matching
%%%%%%%%%%%%%%%%%%%%%%%%%
Mat_diffs = 100*ones(size_omega,size_sigma_a,size_sigma_g);
num_cells = size_omega*size_sigma_a*size_sigma_g;
for cont_omega = 1:size_omega
cont_omega
size_omega
end
for cont_sigma_a = 1:size_sigma_a
for cont_sigma_g = 1:size_sigma_g
OMEGA_a = omega_grid(1,cont_omega);
OMEGA_g = OMEGA_a;
SIGMA_a = sigma_a_grid(1,cont_sigma_a);
SIGMA_g = sigma_g_grid(1,cont_sigma_g);
%%%%%%%%%%%%%%%%%%%
% Solving Model
%%%%%%%%%%%%%%%%%%%
[fx,fxp,fy,fyp] = benchmark_model(RHOA,SIGMA_a,RHOG,SIGMA_g,TAUN,TAUY,E,GAMMA,OMEGA_a,OMEGA_g,DELTA_F,DELTA_I,ALPHA_I,CHIF,CHII);
[gx,hx] = gx_hx(fy,fx,fyp,fxp);
% Analytical-Model-Based HP Second Moments
varshock = zeros(size(hx,1));
varshock(4,4) = SIGMA_a^2;
varshock(5,5) = SIGMA_g^2;
[var_y,var_x] = mom(gx,hx,varshock,0); %Compute Variance/covariance matrix
[var_y1,var_x1] = mom(gx,hx,varshock,1); %First-order autocovariance
table_modelmoms(1:13,1) = (diag(var_y(1:13,1:13)).^(1/2))*100;
table_modelmoms(1:13,2) = diag(var_y1(1:13,1:13)) ./ diag(var_y(1:13,1:13));
table_modelmoms(1:13,3) = var_y(1:13,1) ./ (diag(var_y(1:13,1:13)).^(1/2)) /var_y(1,1)^(1/2);
% Table of serial correlation Labor and Output
corr_yFt_hIt = var_y(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
corr_yFt_hIt_plus1 = var_y1(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y2,var_x2] = mom(gx,hx,varshock,2);
corr_yFt_hIt_plus2 = var_y2(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y3,var_x3] = mom(gx,hx,varshock,3);
corr_yFt_hIt_plus3 = var_y3(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y4,var_x4] = mom(gx,hx,varshock,4);
corr_yFt_hIt_plus4 = var_y4(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym1,var_xm1] = mom(gx,hx,varshock,-1);
corr_yFt_hIt_minus1 = var_ym1(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym2,var_xm2] = mom(gx,hx,varshock,-2);
corr_yFt_hIt_minus2 = var_ym2(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym3,var_xm3] = mom(gx,hx,varshock,-3);
corr_yFt_hIt_minus3 = var_ym3(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym4,var_xm4] = mom(gx,hx,varshock,-4);
corr_yFt_hIt_minus4 = var_ym4(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
MOMS_hp_model = ...
[table_modelmoms(1,1); % std(yF)*100 %1
table_modelmoms(2,1)/table_modelmoms(1,1); % std(cF)/std(yF) %2
table_modelmoms(3,1)/table_modelmoms(1,1); % std(ivF)/std(yF); %3
table_modelmoms(7,1)/table_modelmoms(1,1); % std(tby)/std(yF); %4
table_modelmoms(1,2); % corr(yF_t,yF_t-1);%5
table_modelmoms(2,3); % corr(yF,cF); %6
table_modelmoms(3,3); % corr(yF,ivF); %7
table_modelmoms(7,3); % corr(yF,tby); %8
table_modelmoms(6,1)/table_modelmoms(1,1); % std(h)/std(yF); %9
table_modelmoms(6,3); % corr(yF,h); %10
table_modelmoms(4,1)/table_modelmoms(1,1); % std(hF)/std(yF); %11
table_modelmoms(5,1)/table_modelmoms(1,1); % std(hI)/std(yF); %12
table_modelmoms(4,3); % corr(yF,hF); %13
table_modelmoms(5,3); % corr(yF,hI); %14
corr_yFt_hIt_minus4; % corr(yF_t,hI_t-4);%15
corr_yFt_hIt_minus3; % corr(yF_t,hI_t-3);%16
corr_yFt_hIt_minus2; % corr(yF_t,hI_t-2);%17
corr_yFt_hIt_minus1; % corr(yF_t,hI_t-1);%18
corr_yFt_hIt_plus1; % corr(yF_t,hI_t+1);%19
corr_yFt_hIt_plus2; % corr(yF_t,hI_t+2);%20
corr_yFt_hIt_plus3; % corr(yF_t,hI_t+3);%21
corr_yFt_hIt_plus4]; % corr(yF_t,hI_t+4);%22
diff = ((MOMS_hp_model-MOMS_hp_data_all)./MOMS_hp_data_all);
diff_1 = diff'*diff;
Mat_diffs(cont_omega,cont_sigma_a,cont_sigma_g) = diff_1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computing moments of the lowest value in the grid.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[mini,posi] = min(Mat_diffs(:));
[loc_omega,loc_sigma_a,loc_sigma_g] = ind2sub([size_omega size_sigma_a size_sigma_g],posi);
omega_a_est = omega_grid(loc_omega); % Calibrated omega presented in upper panel of Table 6, Column labelled "extension 2".
omega_g_est = omega_a_est;
sigma_a_est = sigma_a_grid(loc_sigma_a); % Calibrated sigma_a presented in upper panel of Table 6, Column labelled "extension 2".
sigma_g_est = sigma_g_grid(loc_sigma_g); % Calibrated sigma_g presented in upper panel of Table 6, Column labelled "extension 2".
[fx,fxp,fy,fyp] = benchmark_model(RHOA,SIGMA_a,RHOG,SIGMA_g,TAUN,TAUY,E,GAMMA,OMEGA_a,OMEGA_g,DELTA_F,DELTA_I,ALPHA_I,CHIF,CHII);
[gx,hx] = gx_hx(fy,fx,fyp,fxp);
% Analytical-Model-Based HP Second Moments
% Table of second moments
varshock = zeros(size(hx,1));
varshock(4,4) = sigma_a_est^2;
varshock(5,5) = sigma_g_est^2;
[var_y,var_x] = mom(gx,hx,varshock,0); %Compute Variance/covariance matrix
[var_y1,var_x1] = mom(gx,hx,varshock,1); %First-order autocovariance
table_modelmoms(1:13,1) = (diag(var_y(1:13,1:13)).^(1/2))*100;
table_modelmoms(1:13,2) = diag(var_y1(1:13,1:13)) ./ diag(var_y(1:13,1:13));
table_modelmoms(1:13,3) = var_y(1:13,1) ./ (diag(var_y(1:13,1:13)).^(1/2)) /var_y(1,1)^(1/2);
% Table of serial correlation
corr_yFt_hIt = var_y(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
corr_yFt_hIt_plus1 = var_y1(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y2,var_x2] = mom(gx,hx,varshock,2); %Second-order autocovariance
corr_yFt_hIt_plus2 = var_y2(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y3,var_x3] = mom(gx,hx,varshock,3); %Third-order autocovariance
corr_yFt_hIt_plus3 = var_y3(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y4,var_x4] = mom(gx,hx,varshock,4); %Fourth-order autocovariance
corr_yFt_hIt_plus4 = var_y4(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym1,var_xm1] = mom(gx,hx,varshock,-1);
corr_yFt_hIt_minus1 = var_ym1(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym2,var_xm2] = mom(gx,hx,varshock,-2);
corr_yFt_hIt_minus2 = var_ym2(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym3,var_xm3] = mom(gx,hx,varshock,-3);
corr_yFt_hIt_minus3 = var_ym3(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym4,var_xm4] = mom(gx,hx,varshock,-4);
corr_yFt_hIt_minus4 = var_ym4(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
MOMS_hp_model_est = ...
[table_modelmoms(1,1); % std(yF)*100 %1
table_modelmoms(2,1)/table_modelmoms(1,1); % std(cF)/std(yF) %2
table_modelmoms(3,1)/table_modelmoms(1,1); % std(ivF)/std(yF); %3
table_modelmoms(7,1)/table_modelmoms(1,1); % std(tby)/std(yF); %4
table_modelmoms(1,2); % corr(yF_t,yF_t-1);%5
table_modelmoms(2,3); % corr(yF,cF); %6
table_modelmoms(3,3); % corr(yF,ivF); %7
table_modelmoms(7,3); % corr(yF,tby); %8
table_modelmoms(6,1)/table_modelmoms(1,1); % std(h)/std(yF); %9
table_modelmoms(6,3); % corr(yF,h); %10
table_modelmoms(4,1)/table_modelmoms(1,1); % std(hF)/std(yF); %11
table_modelmoms(5,1)/table_modelmoms(1,1); % std(hI)/std(yF); %12
table_modelmoms(4,3); % corr(yF,hF); %13
table_modelmoms(5,3); % corr(yF,hI); %14
corr_yFt_hIt_minus4; % corr(yF_t,hI_t-4);%15
corr_yFt_hIt_minus3; % corr(yF_t,hI_t-3);%16
corr_yFt_hIt_minus2; % corr(yF_t,hI_t-2);%17
corr_yFt_hIt_minus1; % corr(yF_t,hI_t-1);%18
corr_yFt_hIt_plus1; % corr(yF_t,hI_t+1);%19
corr_yFt_hIt_plus2; % corr(yF_t,hI_t+2);%20
corr_yFt_hIt_plus3; % corr(yF_t,hI_t+3);%21
corr_yFt_hIt_plus4]; % corr(yF_t,hI_t+4);%22
MOMS_comparison_2moms = [MOMS_hp_data_all(1:14,1),MOMS_hp_model_est(1:14,1)]; % Final results reported in Table 3 / Column 3
MOMS_comparison_serial = [MOMS_hp_data_corrs_yF_hI,...
[corr_yFt_hIt_minus4;...
corr_yFt_hIt_minus3;...
corr_yFt_hIt_minus2;...
corr_yFt_hIt_minus1;...
corr_yFt_hIt;...
corr_yFt_hIt_plus1;...
corr_yFt_hIt_plus2;...
corr_yFt_hIt_plus3;...
corr_yFt_hIt_plus4;...
]];
% Final results reported in Figure 8. In workspace we
% call this as serial_hI_AG
time_serial = [-4,-3,-2,-1,0,1,2,3,4]';
figure(1)
plot(time_serial,MOMS_comparison_serial(:,1),'--r',time_serial,MOMS_comparison_serial(:,2),'b',...
time_serial,MOMS_comparison_serial(:,3),'--r',time_serial,MOMS_comparison_serial(:,4))
Answers (1)
Star Strider
on 4 Feb 2018
0 votes
What you got was not an error. It was the result of the fsolve optimisation (or parameter estimation or whatever you are doing), and reported it being successful. The 1 returned as ‘exitflag’ means ‘Equation solved. First-order optimality is small.’ The reason you also got a 0 for ‘exitflag’ (number of iterations or function evaluations exceeded the limits set for them) apparently arises from another fsolve call you have not explained.
Use the parameter values it returned with 1 for ‘exitflag’ for whatever you need to do, because they are optimal.
2 Comments
AYSEGUL DILSIZ
on 4 Feb 2018
Star Strider
on 4 Feb 2018
My pleasure.
The are in the output from fsolve. I do not see anywhere in your code where you called fsolve.
If you called it like this:
x = fsolve(fun,x0)
the optimised parameters are the ‘x’ vector.
Categories
Find more on Common Operations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!