diff --git a/doc/dynare.texi b/doc/dynare.texi index 49729fdf9..8644a0db5 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4924,6 +4924,9 @@ Default value is @code{0}. forecast covariance matrices. @item filter_step_ahead = [@var{INTEGER1}:@var{INTEGER2}] +See below. + +@item filter_step_ahead = [@var{INTEGER1} @var{INTEGER2} @dots{}] @anchor{filter_step_ahead} @vindex oo_.FilteredVariablesKStepAhead @vindex oo_.FilteredVariablesKStepAheadVariances @@ -5186,7 +5189,7 @@ Variable set by the @code{estimation} command, if it is used with the @defvr {MATLAB/Octave variable} oo_.FilteredVariablesKStepAhead Variable set by the @code{estimation} command, if it is used with the -@code{filter_step_ahead} option. +@code{filter_step_ahead} option. The k-steps are stored along the rows while the columns indicate the respective variables. The third dimension of the array provides the observation for which the forecast has been made. For example, if @code{filter_step_ahead=[1 2 4]} and @code{nobs=200}, the element (3,5,204) stores the four period ahead filtered value of variable 5 computed at time t=200 for time t=204. The periods at the beginning and end of the sample for which no forecasts can be made, e.g. entries (1,5,1) and (1,5,204) in the example, are set to zero. @end defvr @defvr {MATLAB/Octave variable} oo_.FilteredVariablesKStepAheadVariances @@ -5194,6 +5197,16 @@ Variable set by the @code{estimation} command, if it is used with the @code{filter_step_ahead} option. @end defvr +@defvr {MATLAB/Octave variable} oo_.Filtered_Variables_X_step_ahead +Variable set by the @code{estimation} command, if it is used with the @code{filter_step_ahead} option in the context of Bayesian estimation. Fields are of the form: +@example +@code{oo_.Filtered_Variables_X_step_ahead.@var{VARIABLE_NAME}} +@end example +The nth entry stores the k-step ahead filtered variable computed at time n for time n+k. +@end defvr + + + @defvr {MATLAB/Octave variable} oo_.PosteriorIRF.dsge Variable set by the @code{estimation} command, if it is used with the @code{bayesian_irf} option. Fields are of the form: diff --git a/matlab/pm3.m b/matlab/pm3.m index fa2bb0ec2..567ba71da 100644 --- a/matlab/pm3.m +++ b/matlab/pm3.m @@ -56,29 +56,67 @@ if options_.TeX end Mean = zeros(n2,nvar); Median = zeros(n2,nvar); -Std = zeros(n2,nvar); +Var = zeros(n2,nvar); Distrib = zeros(9,n2,nvar); HPD = zeros(2,n2,nvar); fprintf(['Estimation::mcmc: ' tit1 '\n']); stock1 = zeros(n1,n2,B); k = 0; +filter_step_ahead_indicator=0; for file = 1:ifil load([DirectoryName '/' M_.fname var_type int2str(file)]); if size(size(stock),2) == 4 - stock = squeeze(stock(1,:,1:n2,:)); + filter_step_ahead_indicator=1; + stock_filter_step_ahead=zeros(n1,n2,size(stock,4),length(options_.filter_step_ahead)); + for ii=1:length(options_.filter_step_ahead) + K_step_ahead=options_.filter_step_ahead(ii); + stock_filter_step_ahead(:,:,:,ii)=stock(ii,:,1+K_step_ahead:n2+K_step_ahead,:); + end + stock = squeeze(stock(1,:,1+1:1+n2,:)); %1 step ahead starts at entry 2 end k = k(end)+(1:size(stock,3)); stock1(:,:,k) = stock; + if filter_step_ahead_indicator + stock1_filter_step_ahead(:,:,k,:) = stock_filter_step_ahead; + end end clear stock +if filter_step_ahead_indicator + clear stock_filter_step_ahead + filter_steps=length(options_.filter_step_ahead); + Mean_filter_step_ahead = zeros(filter_steps,nvar,n2); + Median_filter_step_ahead = zeros(filter_steps,nvar,n2); + Var_filter_step_ahead = zeros(filter_steps,nvar,n2); + Distrib_filter_step_ahead = zeros(9,filter_steps,nvar,n2); + HPD_filter_step_ahead = zeros(2,filter_steps,nvar,n2); +end + tmp =zeros(B,1); for i = 1:nvar for j = 1:n2 [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ... posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig); + if filter_step_ahead_indicator + for K_step = 1:length(options_.filter_step_ahead) + [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ... + posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig); + end + end end end clear stock1 +if filter_step_ahead_indicator %write matrices corresponding to ML + clear stock1_filter_step_ahead + FilteredVariablesKStepAhead=zeros(length(options_.filter_step_ahead),nvar,n2+max(options_.filter_step_ahead)); + FilteredVariablesKStepAheadVariances=zeros(length(options_.filter_step_ahead),nvar,n2+max(options_.filter_step_ahead)); + for K_step = 1:length(options_.filter_step_ahead) + FilteredVariablesKStepAhead(K_step,:,1+options_.filter_step_ahead(K_step):n2+options_.filter_step_ahead(K_step))=Mean_filter_step_ahead(K_step,:,:); + FilteredVariablesKStepAheadVariances(K_step,:,1+options_.filter_step_ahead(K_step):n2+options_.filter_step_ahead(K_step))=Mean_filter_step_ahead(K_step,:,:); + end + oo_.FilteredVariablesKStepAhead=FilteredVariablesKStepAhead; + oo_.FilteredVariablesKStepAheadVariances=FilteredVariablesKStepAheadVariances; +end + for i = 1:nvar name = deblank(names1(SelecVariables(i),:)); eval(['oo_.' name3 '.Mean.' name ' = Mean(:,i);']); @@ -87,6 +125,17 @@ for i = 1:nvar eval(['oo_.' name3 '.deciles.' name ' = Distrib(:,:,i);']); eval(['oo_.' name3 '.HPDinf.' name ' = HPD(1,:,i);']); eval(['oo_.' name3 '.HPDsup.' name ' = HPD(2,:,i);']); + if filter_step_ahead_indicator + for K_step = 1:length(options_.filter_step_ahead) + name4=['Filtered_Variables_',num2str(K_step),'_step_ahead']; + eval(['oo_.' name4 '.Mean.' name ' = squeeze(Mean_filter_step_ahead(K_step,i,:));']); + eval(['oo_.' name4 '.Median.' name ' = squeeze(Median_filter_step_ahead(K_step,i,:));']); + eval(['oo_.' name4 '.Var.' name ' = squeeze(Var_filter_step_ahead(K_step,i,:));']); + eval(['oo_.' name4 '.deciles.' name ' = squeeze(Distrib_filter_step_ahead(:,K_step,i,:));']); + eval(['oo_.' name4 '.HPDinf.' name ' = squeeze(HPD_filter_step_ahead(1,K_step,i,:));']); + eval(['oo_.' name4 '.HPDsup.' name ' = squeeze(HPD_filter_step_ahead(2,K_step,i,:));']); + end + end end %% %% Finally I build the plots. diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index f62b73019..31de4cab3 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -300,7 +300,7 @@ if options_.filtered_vars '',varlist,M_.endo_names_tex,M_.endo_names,... varlist,'UpdatedVariables',DirectoryName, ... '_update'); - pm3(endo_nbr,gend+1,ifil(4),B,'One step ahead forecast (filtered variables)',... + pm3(endo_nbr,gend,ifil(4),B,'One step ahead forecast (filtered variables)',... '',varlist,M_.endo_names_tex,M_.endo_names,... varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead'); end diff --git a/tests/filter_step_ahead/fs2000_filter_step_ahead_ML.mod b/tests/filter_step_ahead/fs2000_filter_step_ahead_ML.mod new file mode 100644 index 000000000..c00a80e15 --- /dev/null +++ b/tests/filter_step_ahead/fs2000_filter_step_ahead_ML.mod @@ -0,0 +1,112 @@ +/* + * This file replicates the estimation of the cash in advance model described + * Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models", + * Journal of Applied Econometrics, 15(6), 645-670. + * + * The data are in file "fsdat_simul.m", and have been artificially generated. + * They are therefore different from the original dataset used by Schorfheide. + * + * The equations are taken from J. Nason and T. Cogley (1994): "Testing the + * implications of long-run neutrality for monetary business cycle models", + * Journal of Applied Econometrics, 9, S37-S70. + * Note that there is an initial minus sign missing in equation (A1), p. S63. + * + * This implementation was written by Michel Juillard. Please note that the + * following copyright notice only applies to this Dynare implementation of the + * model. + */ + +/* + * Copyright (C) 2004-2010 Dynare Team + * + * This file is part of Dynare. + * + * Dynare 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 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 . + */ + +var m P c e W R k d n l gy_obs gp_obs y dA; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +initval; +k = 6; +m = mst; +P = 2.25; +c = 0.45; +e = 1; +W = 4; +R = 1.02; +d = 0.85; +n = 0.19; +l = 0.86; +y = 0.6; +gy_obs = exp(gam); +gp_obs = exp(-gam); +dA = exp(gam); +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +check; + +estimated_params; +alp, 0.356; +bet, 0.993; +gam, 0.0085; +end; + +varobs gp_obs gy_obs; + +estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, filter_step_ahead = [1 4 8 12], forecast=20,smoother,filtered_vars) m P c; + + +/* + * The following lines were used to generate the data file. If you want to + * generate another random data file, comment the "estimation" line and uncomment + * the following lines. + */ + +//stoch_simul(periods=200, order=1); +//datatomfile('fsdat_simul', char('gy_obs', 'gp_obs')); diff --git a/tests/filter_step_ahead/fs2000_filter_step_ahead_ML_steadystate.m b/tests/filter_step_ahead/fs2000_filter_step_ahead_ML_steadystate.m new file mode 100644 index 000000000..2ac140da5 --- /dev/null +++ b/tests/filter_step_ahead/fs2000_filter_step_ahead_ML_steadystate.m @@ -0,0 +1,73 @@ +% computes the steady state of fs2000 analyticaly +% largely inspired by the program of F. Schorfheide + +% Copyright (C) 2004-2010 Dynare Team +% +% This file is part of Dynare. +% +% Dynare 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 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 . + +function [ys,check] = fs2000_steadystate(ys,exe) + global M_ + + alp = M_.params(1); + bet = M_.params(2); + gam = M_.params(3); + mst = M_.params(4); + rho = M_.params(5); + psi = M_.params(6); + del = M_.params(7); + + check = 0; + + dA = exp(gam); + gst = 1/dA; + m = mst; + + khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); + xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); + nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); + n = xist/(nust+xist); + P = xist + nust; + k = khst*n; + + l = psi*mst*n/( (1-psi)*(1-n) ); + c = mst/P; + d = l - mst + 1; + y = k^alp*n^(1-alp)*gst^alp; + R = mst/bet; + W = l/n; + ist = y-c; + q = 1 - d; + + e = 1; + + gp_obs = m/dA; + gy_obs = dA; + + ys =[ +m +P +c +e +W +R +k +d +n +l +gy_obs +gp_obs +y +dA ]; diff --git a/tests/filter_step_ahead/fs2000_filter_step_ahead_bayesian.mod b/tests/filter_step_ahead/fs2000_filter_step_ahead_bayesian.mod new file mode 100644 index 000000000..928029f64 --- /dev/null +++ b/tests/filter_step_ahead/fs2000_filter_step_ahead_bayesian.mod @@ -0,0 +1,118 @@ +/* + * This file replicates the estimation of the cash in advance model described + * Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models", + * Journal of Applied Econometrics, 15(6), 645-670. + * + * The data are in file "fsdat_simul.m", and have been artificially generated. + * They are therefore different from the original dataset used by Schorfheide. + * + * The equations are taken from J. Nason and T. Cogley (1994): "Testing the + * implications of long-run neutrality for monetary business cycle models", + * Journal of Applied Econometrics, 9, S37-S70. + * Note that there is an initial minus sign missing in equation (A1), p. S63. + * + * This implementation was written by Michel Juillard. Please note that the + * following copyright notice only applies to this Dynare implementation of the + * model. + */ + +/* + * Copyright (C) 2004-2010 Dynare Team + * + * This file is part of Dynare. + * + * Dynare 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 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 . + */ + +var m P c e W R k d n l gy_obs gp_obs y dA; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +initval; +k = 6; +m = mst; +P = 2.25; +c = 0.45; +e = 1; +W = 4; +R = 1.02; +d = 0.85; +n = 0.19; +l = 0.86; +y = 0.6; +gy_obs = exp(gam); +gp_obs = exp(-gam); +dA = exp(gam); +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +check; + +estimated_params; +alp, beta_pdf, 0.356, 0.02; +bet, beta_pdf, 0.993, 0.002; +gam, normal_pdf, 0.0085, 0.003; +mst, normal_pdf, 1.0002, 0.007; +rho, beta_pdf, 0.129, 0.223; +psi, beta_pdf, 0.65, 0.05; +del, beta_pdf, 0.01, 0.005; +stderr e_a, inv_gamma_pdf, 0.035449, inf; +stderr e_m, inv_gamma_pdf, 0.008862, inf; +end; + +varobs gp_obs gy_obs; + +estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=2000, mh_nblocks=1, mh_jscale=0.8,filter_step_ahead = [1 4 8 12], forecast=20,smoother,filtered_vars) m P c; + + +/* + * The following lines were used to generate the data file. If you want to + * generate another random data file, comment the "estimation" line and uncomment + * the following lines. + */ + +//stoch_simul(periods=200, order=1); +//datatomfile('fsdat_simul', char('gy_obs', 'gp_obs')); diff --git a/tests/filter_step_ahead/fs2000_filter_step_ahead_bayesian_steadystate.m b/tests/filter_step_ahead/fs2000_filter_step_ahead_bayesian_steadystate.m new file mode 100644 index 000000000..2ac140da5 --- /dev/null +++ b/tests/filter_step_ahead/fs2000_filter_step_ahead_bayesian_steadystate.m @@ -0,0 +1,73 @@ +% computes the steady state of fs2000 analyticaly +% largely inspired by the program of F. Schorfheide + +% Copyright (C) 2004-2010 Dynare Team +% +% This file is part of Dynare. +% +% Dynare 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 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 . + +function [ys,check] = fs2000_steadystate(ys,exe) + global M_ + + alp = M_.params(1); + bet = M_.params(2); + gam = M_.params(3); + mst = M_.params(4); + rho = M_.params(5); + psi = M_.params(6); + del = M_.params(7); + + check = 0; + + dA = exp(gam); + gst = 1/dA; + m = mst; + + khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); + xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); + nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); + n = xist/(nust+xist); + P = xist + nust; + k = khst*n; + + l = psi*mst*n/( (1-psi)*(1-n) ); + c = mst/P; + d = l - mst + 1; + y = k^alp*n^(1-alp)*gst^alp; + R = mst/bet; + W = l/n; + ist = y-c; + q = 1 - d; + + e = 1; + + gp_obs = m/dA; + gy_obs = dA; + + ys =[ +m +P +c +e +W +R +k +d +n +l +gy_obs +gp_obs +y +dA ]; diff --git a/tests/filter_step_ahead/fsdat_simul.m b/tests/filter_step_ahead/fsdat_simul.m new file mode 100644 index 000000000..56c0e4cd5 --- /dev/null +++ b/tests/filter_step_ahead/fsdat_simul.m @@ -0,0 +1,416 @@ +% Generated data, used by fs2000.mod + +gy_obs =[ + 1.0030045 + 1.0002599 + 0.99104664 + 1.0321162 + 1.0223545 + 1.0043614 + 0.98626929 + 1.0092127 + 1.0357197 + 1.0150827 + 1.0051548 + 0.98465775 + 0.99132132 + 0.99904153 + 1.0044641 + 1.0179198 + 1.0113462 + 0.99409421 + 0.99904293 + 1.0448336 + 0.99932433 + 1.0057004 + 0.99619787 + 1.0267504 + 1.0077645 + 1.0058026 + 1.0025891 + 0.9939097 + 0.99604693 + 0.99908569 + 1.0151094 + 0.99348134 + 1.0039124 + 1.0145805 + 0.99800868 + 0.98578138 + 1.0065771 + 0.99843919 + 0.97979062 + 0.98413351 + 0.96468174 + 1.0273857 + 1.0225211 + 0.99958667 + 1.0111157 + 1.0099585 + 0.99480311 + 1.0079265 + 0.98924573 + 1.0070613 + 1.0075706 + 0.9937151 + 1.0224711 + 1.0018891 + 0.99051863 + 1.0042944 + 1.0184055 + 0.99419508 + 0.99756624 + 1.0015983 + 0.9845772 + 1.0004407 + 1.0116237 + 0.9861885 + 1.0073094 + 0.99273355 + 1.0013224 + 0.99777979 + 1.0301686 + 0.96809556 + 0.99917088 + 0.99949253 + 0.96590004 + 1.0083938 + 0.96662298 + 1.0221454 + 1.0069792 + 1.0343996 + 1.0066531 + 1.0072525 + 0.99743563 + 0.99723703 + 1.000372 + 0.99013917 + 1.0095223 + 0.98864268 + 0.98092242 + 0.98886488 + 1.0030341 + 1.01894 + 0.99155059 + 0.99533235 + 0.99734316 + 1.0047356 + 1.0082737 + 0.98425116 + 0.99949212 + 1.0055899 + 1.0065075 + 0.99385069 + 0.98867975 + 0.99804843 + 1.0184038 + 0.99301902 + 1.0177222 + 1.0051924 + 1.0187852 + 1.0098985 + 1.0097172 + 1.0145811 + 0.98721038 + 1.0361722 + 1.0105821 + 0.99469309 + 0.98626785 + 1.013871 + 0.99858924 + 0.99302637 + 1.0042186 + 0.99623745 + 0.98545708 + 1.0225435 + 1.0011861 + 1.0130321 + 0.97861347 + 1.0228193 + 0.99627435 + 1.0272779 + 1.0075172 + 1.0096762 + 1.0129306 + 0.99966549 + 1.0262882 + 1.0026914 + 1.0061475 + 1.009523 + 1.0036127 + 0.99762992 + 0.99092634 + 1.0058469 + 0.99887292 + 1.0060653 + 0.98673557 + 0.98895709 + 0.99111967 + 0.990118 + 0.99788054 + 0.97054709 + 1.0099157 + 1.0107431 + 0.99518695 + 1.0114048 + 0.99376019 + 1.0023369 + 0.98783327 + 1.0051727 + 1.0100462 + 0.98607387 + 1.0000064 + 0.99692442 + 1.012225 + 0.99574078 + 0.98642833 + 0.99008207 + 1.0197359 + 1.0112849 + 0.98711069 + 0.99402748 + 1.0242141 + 1.0135349 + 0.99842505 + 1.0130714 + 0.99887044 + 1.0059058 + 1.0185998 + 1.0073314 + 0.98687706 + 1.0084551 + 0.97698964 + 0.99482714 + 1.0015302 + 1.0105331 + 1.0261767 + 1.0232822 + 1.0084176 + 0.99785167 + 0.99619733 + 1.0055223 + 1.0076326 + 0.99205461 + 1.0030587 + 1.0137012 + 1.0145878 + 1.0190297 + 1.0000681 + 1.0153894 + 1.0140649 + 1.0007236 + 0.97961463 + 1.0125257 + 1.0169503 + 1.0197363 + 1.0221185 + +]; + +gp_obs =[ + 1.0079715 + 1.0115853 + 1.0167502 + 1.0068957 + 1.0138189 + 1.0258364 + 1.0243817 + 1.017373 + 1.0020171 + 1.0003742 + 1.0008974 + 1.0104804 + 1.0116393 + 1.0114294 + 0.99932124 + 0.99461459 + 1.0170349 + 1.0051446 + 1.020639 + 1.0051964 + 1.0093042 + 1.007068 + 1.01086 + 0.99590086 + 1.0014883 + 1.0117332 + 0.9990095 + 1.0108284 + 1.0103672 + 1.0036722 + 1.0005124 + 1.0190331 + 1.0130978 + 1.007842 + 1.0285436 + 1.0322054 + 1.0213403 + 1.0246486 + 1.0419306 + 1.0258867 + 1.0156316 + 0.99818589 + 0.9894107 + 1.0127584 + 1.0146882 + 1.0136529 + 1.0340107 + 1.0343652 + 1.02971 + 1.0077932 + 1.0198114 + 1.013971 + 1.0061083 + 1.0089573 + 1.0037926 + 1.0082071 + 0.99498155 + 0.99735772 + 0.98765026 + 1.006465 + 1.0196088 + 1.0053233 + 1.0119974 + 1.0188066 + 1.0029302 + 1.0183459 + 1.0034218 + 1.0158799 + 0.98824798 + 1.0274357 + 1.0168832 + 1.0180641 + 1.0294657 + 0.98864091 + 1.0358326 + 0.99889969 + 1.0178322 + 0.99813566 + 1.0073549 + 1.0215985 + 1.0084245 + 1.0080939 + 1.0157021 + 1.0075815 + 1.0032633 + 1.0117871 + 1.0209276 + 1.0077569 + 0.99680958 + 1.0120266 + 1.0017625 + 1.0138811 + 1.0198358 + 1.0059629 + 1.0115416 + 1.0319473 + 1.0167074 + 1.0116111 + 1.0048627 + 1.0217622 + 1.0125221 + 1.0142045 + 0.99792469 + 0.99823971 + 0.99561547 + 0.99850373 + 0.9898464 + 1.0030963 + 1.0051373 + 1.0004213 + 1.0144117 + 0.97185592 + 0.9959518 + 1.0073529 + 1.0051603 + 0.98642572 + 0.99433423 + 1.0112131 + 1.0007695 + 1.0176867 + 1.0134363 + 0.99926191 + 0.99879835 + 0.99878754 + 1.0331374 + 1.0077797 + 1.0127221 + 1.0047393 + 1.0074106 + 0.99784213 + 1.0056495 + 1.0057708 + 0.98817494 + 0.98742176 + 0.99930555 + 1.0000687 + 1.0129754 + 1.009529 + 1.0226731 + 1.0149534 + 1.0164295 + 1.0239469 + 1.0293458 + 1.026199 + 1.0197525 + 1.0126818 + 1.0054473 + 1.0254423 + 1.0069461 + 1.0153135 + 1.0337515 + 1.0178187 + 1.0240469 + 1.0079489 + 1.0186953 + 1.0008628 + 1.0113799 + 1.0140118 + 1.0168007 + 1.011441 + 0.98422774 + 0.98909729 + 1.0157859 + 1.0151586 + 0.99756232 + 0.99497777 + 1.0102841 + 1.0221659 + 0.9937759 + 0.99877193 + 1.0079433 + 0.99667692 + 1.0095959 + 1.0128804 + 1.0156949 + 1.0111951 + 1.0228887 + 1.0122083 + 1.0190197 + 1.0074927 + 1.0268096 + 0.99689352 + 0.98948474 + 1.0024938 + 1.0105543 + 1.014116 + 1.0141217 + 1.0056504 + 1.0101026 + 1.0105069 + 0.99619053 + 1.0059439 + 0.99449473 + 0.99482458 + 1.0037702 + 1.0068087 + 0.99575975 + 1.0030815 + 1.0334014 + 0.99879386 + 0.99625634 + 1.0171195 + 0.99233844 + +]; +