diff --git a/doc/dynare.texi b/doc/dynare.texi index 7265b4aee..d31f2be57 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -9733,7 +9733,7 @@ Instantiates an empty @dseries object, with, if defined, an initial date given b @sp 1 -@deftypefn {dseries} dseries (@var{FILENAME}) +@deftypefn {dseries} dseries (@var{FILENAME}[, @var{INITIAL_DATE}]) Instantiates and populates a @dseries object with a data file specified by @var{FILENAME}, a string passed as input. Valid file types are @file{.m} file, @file{.mat} file, @file{.csv} file, and @file{.xls} file. A typical @file{.m} file will have the following form: @@ -9746,7 +9746,7 @@ azert = randn(100,1); yuiop = randn(100,1); @end example -If a @file{.mat} file is used instead, it should provide the same informations. Note that the @code{INIT__} variable can be either a @dates object or a string which could be used to instantiate the same @dates object. +If a @file{.mat} file is used instead, it should provide the same informations. Note that the @code{INIT__} variable can be either a @dates object or a string which could be used to instantiate the same @dates object. If @code{INIT__} is not provided in the @file{.mat} or @file{.m} file, the initial is by default set equal to @code{dates('1Y')}. If a second input argument is passed to the constructor, @dates object @var{INITIAL_DATE}, the initial date defined in @var{FILENAME} is reset to @var{INITIAL_DATE}. This is typically usefull if @code{INIT__} is not provided in the data file. @end deftypefn diff --git a/matlab/@dates/dates.m b/matlab/@dates/dates.m index 99093bfe7..c1c80fea2 100644 --- a/matlab/@dates/dates.m +++ b/matlab/@dates/dates.m @@ -45,7 +45,7 @@ function dd = dates(varargin) % --*-- Unitary tests --*-- %! @end deftypefn %@eod: -% Copyright (C) 2011-2013 Dynare Team +% Copyright (C) 2011-2014 Dynare Team % % This file is part of Dynare. % @@ -62,6 +62,7 @@ function dd = dates(varargin) % --*-- Unitary tests --*-- % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +% Initialization. if nargin>0 && ischar(varargin{1}) && isequal(varargin{1},'initialize') dd = struct('ndat', 0, 'freq', NaN(0), 'time', NaN(0,2)); dd = class(dd,'dates'); @@ -71,61 +72,98 @@ end dd = evalin('base','emptydatesobject'); -switch nargin - case 0 - % Returns an empty object +if isequal(nargin, 0) + % Return an empty dates obect return - case 1 - if isa(varargin{1},'dates') - % Returns a copy of the input argument - dd = varargin{1}; - elseif isdate(varargin{1}) - date = string2date(varargin{1}); - dd.ndat = 1; - dd.freq = date.freq; - dd.time = date.time; - elseif isfreq(varargin{1}) - % Instantiate an empty dates object (only set frequency) - if ischar(varargin{1}) - dd.freq = string2freq(varargin{1}); - else - dd.freq = varargin{1}; - end +end + +if all(cellfun(@isdates, varargin)) + % Concatenates dates in a dates object. + dd = horzcat(varargin{:}); + return +end + +if all(cellfun(@isstringdate,varargin)) + % Concatenates dates in a dates object. + tmp = cellfun(@string2date,varargin); + if all([tmp.freq]-tmp(1).freq==0) + dd.freq = tmp(1).freq; else - error('dates:: Wrong calling sequence of the constructor!') + error('dates::dates: Wrong calling sequence of the constructor! All dates must have common frequency.') end - otherwise - if isdate(varargin{1}) - dd.ndat = nargin; - dd.time = NaN(dd.ndat,2); - date = string2date(varargin{1}); - dd.freq = date.freq; - dd.time(1,:) = date.time; - elseif isdates(varargin{1}) - dd = horzcat(varargin{:}); - return - elseif isfreq(varargin{1}) - S.type = '()'; - S.subs = varargin; - dd = subsref(dd,S); - return + dd.ndat = length(tmp); + dd.time = transpose(reshape([tmp.time],2,dd.ndat)); + return +end + +if isequal(nargin,1) && isfreq(varargin{1}) + % Instantiate an empty dates object (only set frequency) + if ischar(varargin{1}) + dd.freq = string2freq(varargin{1}); else - error('dates::dates: Wrong calling sequence!') + dd.freq = varargin{1}; end - for i=2:dd.ndat - if isdate(varargin{i}) - date = string2date(varargin{i}); - if isequal(date.freq,dd.freq) - dd.time(i,:) = date.time; + return +end + +if isequal(nargin,3) && isfreq(varargin{1}) + if ischar(varargin{1}) + dd.freq = string2freq(varargin{1}); + else + dd.freq = varargin{1}; + end + if (isnumeric(varargin{2}) && isvector(varargin{2}) && isint(varargin{2})) + if isnumeric(varargin{3}) && isvector(varargin{3}) && isint(varargin{3}) + if all(varargin{3}>=1) && all(varargin{3}<=dd.freq) + dd.time = [varargin{2}(:), varargin{3}(:)]; + dd.ndat = size(dd.time,1); else - error(sprintf('dates::dates: Check that all the inputs have the same frequency (see input number %i)!',i)) + error(sprintf('dates::dates: Wrong calling sequence of the constructor! Third input must contain integers between 1 and %i.',dd.freq)) end else - error(sprintf('dates::dates: Input %i has to be a string date!',i)) + error('dates::dates: Wrong calling sequence of the constructor! Third input must be a vector of integers.') + end + else + error('dates::dates: Wrong calling sequence of the constructor! Second input must be a vector of integers.') + end + return +end + +if isequal(nargin,2) && isfreq(varargin{1}) + if ischar(varargin{1}) + dd.freq = string2freq(varargin{1}); + else + dd.freq = varargin{1}; + end + if isequal(dd.freq, 1) + if (isnumeric(varargin{2}) && isvector(varargin{2}) && isint(varargin{2})) + dd.time = [varargin{2}, ones(length(varargin{2}),1)]; + dd.ndat = size(dd.time,1); + return + else + error('dates::dates: Wrong calling sequence of the constructor! Second input must be a vector of integers.') + end + else + if isequal(size(varargin{2},2), 2) + if all(isint(varargin{2}(:,1))) && all(isint(varargin{2}(:,1))) + if all(varargin{2}(:,2)>=1) && all(varargin{2}(:,2)<=dd.freq) + dd.time = [varargin{2}(:,1), varargin{2}(:,2)]; + dd.ndat = size(dd.time,1); + else + error(sprintf('dates::dates: Wrong calling sequence of the constructor! Second column of the last input must contain integers between 1 and %i.',dd.freq)) + end + else + error('dates::dates: Wrong calling sequence! Second input argument must be an array of integers.') + end + else + error('dates::dates: Wrong calling sequence!') end end + return end +error('dates::dates: Wrong calling sequence!') + %@test:1 %$ % Define some dates %$ B1 = '1945Q3'; diff --git a/matlab/@dseries/dseries.m b/matlab/@dseries/dseries.m index 6bcfd2820..502445275 100644 --- a/matlab/@dseries/dseries.m +++ b/matlab/@dseries/dseries.m @@ -10,7 +10,7 @@ function ts = dseries(varargin) % --*-- Unitary tests --*-- %! @sp 2 %! If @code{nargin==0} then an empty dseries object is created. The object can be populated with data subsequently using the overloaded subsref method. %! @sp 2 -%! If @code{nargin==1} and if the input argument is a @ref{dynDate} object, then a dseries object without data is created. This object can be populated with the overload subsref method. +%! If @code{nargin==1} and if the input argument is a @ref{dates} object, then a dseries object without data is created. This object can be populated with the overload subsref method. %! @sp 2 %! If @code{nargin==1} and if the input argument is a string for the name of a csv, m or mat file containing data, then a dseries object is created from these data. %! @sp 2 @@ -54,7 +54,7 @@ function ts = dseries(varargin) % --*-- Unitary tests --*-- %! frequency is unspecified. @var{freq} is equal to 4 if data are on a quaterly basis. @var{freq} is equal to %! 12 if data are on a monthly basis. @var{freq} is equal to 52 if data are on a weekly basis. %! @item init -%! @ref{dynDate} object, initial date of the dataset. +%! @ref{dates} object, initial date of the dataset. %! @end table %! @end deftypefn %@eod: @@ -153,6 +153,12 @@ switch nargin ts.dates = dates(1,1):dates(1,1)+(ts.nobs-1); end case {2,3,4} + if isequal(nargin,2) && ischar(varargin{1}) && isdates(varargin{2}) + % Instantiate dseries object with a data file and force the initial date to be as given by the second input argument. + ds = dseries(varargin{1}); + ts = dseries(ds.data, varargin{2}, ds.name, ds.tex); + return + end a = varargin{1}; b = varargin{2}; if nargin<4 diff --git a/matlab/@dseries/ydiff.m b/matlab/@dseries/ydiff.m index 49b94dc6f..beee58c8a 100644 --- a/matlab/@dseries/ydiff.m +++ b/matlab/@dseries/ydiff.m @@ -151,7 +151,7 @@ end %$ %$ try %$ data = transpose(1:100); -%$ ts = dseries(data,1950,{'A1'},{'A_1'}); +%$ ts = dseries(data,'1950Y',{'A1'},{'A_1'}); %$ ts = ts.ydiff; %$ t(1) = 1; %$ catch diff --git a/matlab/@dynTimeIndex/subsasgn.m b/matlab/@dynTimeIndex/subsasgn.m index b03987c86..101468062 100644 --- a/matlab/@dynTimeIndex/subsasgn.m +++ b/matlab/@dynTimeIndex/subsasgn.m @@ -17,4 +17,4 @@ function val = subsasgn(val, idx, rhs) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -error('dynTimeIndex::subsasgn: Members of dynDate class are private!') \ No newline at end of file +error('dynTimeIndex::subsasgn: Members of dynTimeIndex class are private!') \ No newline at end of file diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index e8e76f97e..d66e77a69 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -59,7 +59,7 @@ R = []; P = []; PK = []; decomp = []; -nobs = size(options_.varobs,1); +nobs = length(options_.varobs); smpl = size(Y,2); M_ = set_all_parameters(xparam1,estim_params_,M_); diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m index 437f8008a..f92c2629e 100644 --- a/matlab/GetPosteriorParametersStatistics.m +++ b/matlab/GetPosteriorParametersStatistics.m @@ -180,17 +180,17 @@ if nvn Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ... posterior_moments(Draws,1,options_.mh_conf_sig); - name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:)); + name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)}; oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density); else try - name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:)); + name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)}; [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type); catch Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ... posterior_moments(Draws,1,options_.mh_conf_sig); - name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:)); + name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)}; oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density); end end diff --git a/matlab/PlotPosteriorDistributions.m b/matlab/PlotPosteriorDistributions.m index 5a622272b..4321103b8 100644 --- a/matlab/PlotPosteriorDistributions.m +++ b/matlab/PlotPosteriorDistributions.m @@ -89,7 +89,7 @@ for i=1:npar eval(['pmod = oo_.posterior_mode.shocks_std.' name ';']) end elseif i <= nvx+nvn - name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i-nvx,1),:)); + name = options_.varobs{estim_params_.nvn_observable_correspondence(i-nvx,1)}; eval(['x1 = oo_.posterior_density.measurement_errors_std.' name '(:,1);']) eval(['f1 = oo_.posterior_density.measurement_errors_std.' name '(:,2);']) eval(['oo_.prior_density.mearsurement_errors_std.' name '(:,1) = x2;']) diff --git a/matlab/PosteriorFilterSmootherAndForecast.m b/matlab/PosteriorFilterSmootherAndForecast.m index 0b3bd1e35..b63656563 100644 --- a/matlab/PosteriorFilterSmootherAndForecast.m +++ b/matlab/PosteriorFilterSmootherAndForecast.m @@ -48,7 +48,7 @@ MaxNumberOfPlotPerFigure = 4;% The square root must be an integer! MaxNumberOfBytes=options_.MaxNumberOfBytes; endo_nbr=M_.endo_nbr; exo_nbr=M_.exo_nbr; -nvobs = size(options_.varobs,1); +nvobs = length(options_.varobs); nn = sqrt(MaxNumberOfPlotPerFigure); iendo = 1:endo_nbr; i_last_obs = gend+(1-M_.maximum_endo_lag:0); @@ -74,9 +74,9 @@ B = 200; MAX_nruns = min(B,ceil(options_.MaxNumberOfBytes/(npar+2)/8)); MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8)); MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8)); -MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8)); +MAX_nerro = min(B,ceil(MaxNumberOfBytes/(length(options_.varobs)*gend)/8)); if naK - MAX_naK = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)* ... + MAX_naK = min(B,ceil(MaxNumberOfBytes/(length(options_.varobs)* ... length(options_.filter_step_ahead)*gend)/8)); end if horizon @@ -272,4 +272,4 @@ if options_.forecast pm3(endo_nbr,horizon+maxlag,ifil6,B,'Forecasted variables (total)',... M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,... 'names2','smooth',MetropolisFolder,'_forc_total') -end \ No newline at end of file +end diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index cf9737f0f..fc1c3bcf5 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -1,11 +1,11 @@ function PosteriorIRF(type) -% Builds posterior IRFs after the MH algorithm. -% -% INPUTS +% Builds posterior IRFs after the MH algorithm. +% +% INPUTS % o type [char] string specifying the joint density of the -% deep parameters ('prior','posterior'). -% -% OUTPUTS +% deep parameters ('prior','posterior'). +% +% OUTPUTS % None (oo_ and plots). % % SPECIAL REQUIREMENTS @@ -34,15 +34,16 @@ function PosteriorIRF(type) % along with Dynare. If not, see . -global options_ estim_params_ oo_ M_ bayestopt_ dataset_ +global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info + % Set the number of periods -if isempty(options_.irf) || ~options_.irf +if isempty(options_.irf) || ~options_.irf options_.irf = 40; end % Set varlist if necessary varlist = options_.varlist; if isempty(varlist) - varlist = options_.varobs; + varlist = char(options_.varobs); end options_.varlist = varlist; nvar = size(varlist,1); @@ -58,7 +59,7 @@ end % Get index of shocks for requested IRFs irf_shocks_indx = getIrfShocksIndx(); - + % Set various parameters & Check or create directories nvx = estim_params_.nvx; nvn = estim_params_.nvn; @@ -68,8 +69,8 @@ np = estim_params_.np ; npar = nvx+nvn+ncx+ncn+np; offset = npar-np; clear('nvx','nvn','ncx','ncn','np'); -nvobs = dataset_.info.nvobs; -gend = dataset_.info.ntobs; +nvobs = dataset_.vobs; +gend = dataset_.nobs; MaxNumberOfPlotPerFigure = 9; nn = sqrt(MaxNumberOfPlotPerFigure); MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8); @@ -111,14 +112,14 @@ else% type = 'prior' B = options_.prior_draws; options_.B = B; end -try +try delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat']) -catch +catch disp('No _IRFs (dsge) files to be deleted!') end -try +try delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat']) -catch +catch disp('No _IRFs (bvar-dsge) files to be deleted!') end irun = 0; @@ -144,10 +145,6 @@ if MAX_nirfs_dsgevar else stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,B); end - [mYY,mXY,mYX,mXX,Ydata,Xdata] = ... - var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,... - options_.dsge_varlag,-1,options_.datafile,options_.varobs,... - options_.xls_sheet,options_.xls_range); NumberOfLags = options_.dsge_varlag; NumberOfLagsTimesNvobs = NumberOfLags*nvobs; if options_.noconstant @@ -174,7 +171,7 @@ localVars.npar = npar; localVars.type=type; if strcmpi(type,'posterior') - while b options_.impulse_responses.plot_threshold - + if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold + name = deblank(varlist(j,:)); texname = deblank(varlist_TeX(j,:)); @@ -413,7 +411,7 @@ if options_.TeX TEXNAMES = char(TEXNAMES,['$' texname '$']); end end - + end fprintf(fidTeX,'\\begin{figure}[H]\n'); for jj = 1:size(TEXNAMES,1) @@ -430,10 +428,10 @@ if options_.TeX fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,' \n'); end - + fprintf(fidTeX,'%% End of TeX file.\n'); fclose(fidTeX); - + end % The others file format are generated in parallel by PosteriorIRF_core2! @@ -453,7 +451,7 @@ if ~isoctave else globalVars = struct('M_',M_, ... 'options_', options_); - + [fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info); end end @@ -463,4 +461,4 @@ end % END parallel code! -fprintf('Estimation::mcmc: Posterior IRFs, done!\n'); \ No newline at end of file +fprintf('Estimation::mcmc: Posterior IRFs, done!\n'); diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m index 1eff387ad..0d0a1a752 100644 --- a/matlab/PosteriorIRF_core1.m +++ b/matlab/PosteriorIRF_core1.m @@ -41,7 +41,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab) % along with Dynare. If not, see . -global options_ estim_params_ oo_ M_ bayestopt_ dataset_ +global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info if nargin<4, whoiam=0; @@ -189,24 +189,24 @@ while fpar1.000000001); end % Get the mean - mu = zeros(1,dataset_.info.nvobs); + mu = zeros(1,dataset_.vobs); % Get rotation if dsge_prior_weight > 0 Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e); @@ -216,24 +216,24 @@ while fpar 1e-10 tmp_dsgevar(:,j) = (irfs(:,j)); end end if IRUN < MAX_nirfs_dsgevar - stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr); + stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr); else - stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr); + stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr); instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ... int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];, eval(['save ' instr]); @@ -242,7 +242,7 @@ while fpar. -global options_ M_ +global options_ M_ if nargin<4, whoiam=0; @@ -87,7 +87,7 @@ OutputFileName={}; subplotnum = 0; for i=fpar:npar, figunumber = 0; - + for j=1:nvar if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold subplotnum = subplotnum+1; @@ -96,7 +96,7 @@ for i=fpar:npar, elseif subplotnum == 1 && ~options_.relative_irf hh = dyn_figure(options_,'Name',['Orthogonalized shock to ' tit(i,:)]); end - + set(0,'CurrentFigure',hh) subplot(nn,nn,subplotnum); if ~MAX_nirfs_dsgevar @@ -108,12 +108,12 @@ for i=fpar:npar, set(h2,'FaceColor',[1 1 1]); set(h2,'BaseValue',min(HPDIRF(:,1,j,i))); plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3) - % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); + % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); box on axis tight xlim([1 options_.irf]); hold off - else + else h1 = area(1:options_.irf,HPDIRF(:,2,j,i)); set(h1,'FaceColor',[.9 .9 .9]); set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))])); @@ -136,9 +136,9 @@ for i=fpar:npar, else if options_.debug fprintf('POSTERIOR_IRF: The IRF of %s to %s is smaller than the irf_plot_threshold of %4.3f and will not be displayed.\n',deblank(varlist(j,:)),tit(i,:),options_.impulse_responses.plot_threshold) - end + end end - + if subplotnum == MaxNumberOfPlotPerFigure || (j == nvar && subplotnum> 0) figunumber = figunumber+1; dyn_saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)],options_); @@ -154,9 +154,8 @@ for i=fpar:npar, % fMessageStatus((i-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab)); dyn_waitbar((i-fpar+1)/(npar-fpar+1),[],waitbarString); end -end% loop over exo_var +end% loop over exo_var myoutput.OutputFileName = OutputFileName; - diff --git a/matlab/ReshapeMatFiles.m b/matlab/ReshapeMatFiles.m index 06acd7372..1ad90fce1 100644 --- a/matlab/ReshapeMatFiles.m +++ b/matlab/ReshapeMatFiles.m @@ -70,7 +70,7 @@ switch type TYPEarray = 4; case 'irf_bvardsge' CAPtype = 'IRF_BVARDSGE'; - TYPEsize = [ options_.irf , size(options_.varobs,1) , M_.exo_nbr ]; + TYPEsize = [ options_.irf , length(options_.varobs) , M_.exo_nbr ]; TYPEarray = 4; case 'smooth' CAPtype = 'SMOOTH'; @@ -82,7 +82,7 @@ switch type TYPEarray = 3; case 'error' CAPtype = 'ERROR'; - TYPEsize = [ size(options_.varobs,1) , options_.nobs ]; + TYPEsize = [ length(options_.varobs) , options_.nobs ]; TYPEarray = 3; case 'innov' CAPtype = 'INNOV'; diff --git a/matlab/bvar_forecast.m b/matlab/bvar_forecast.m index 03100c053..1b88e0ae6 100644 --- a/matlab/bvar_forecast.m +++ b/matlab/bvar_forecast.m @@ -122,7 +122,7 @@ for i = 1:ny dyn_graph=dynare_graph(dyn_graph,[ sims_no_shock_median(:, i) ... sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ... sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ... - options_.varobs(i, :)); + options_.varobs{i}); end dyn_saveas(dyn_graph.fh,[OutputDirectoryName '/' M_.fname '_BVAR_forecast_',num2str(nlags)],options_) @@ -146,8 +146,8 @@ if ~isempty(forecast_data.realized_val) fprintf('RMSE of BVAR(%d):\n', nlags); - for i = 1:size(options_.varobs, 1) - fprintf('%s: %10.4f\n', options_.varobs(i, :), rmse(i)); + for i = 1:length(options_.varobs) + fprintf('%s: %10.4f\n', options_.varobs{i}, rmse(i)); end end @@ -162,8 +162,8 @@ if ~isdir(DirectoryName) end save([ DirectoryName '/simulations.mat'], 'sims_no_shock', 'sims_with_shocks'); -for i = 1:size(options_.varobs, 1) - name = options_.varobs(i, :); +for i = 1:length(options_.varobs) + name = options_.varobs{i}; sims = squeeze(sims_with_shocks(:,i,:)); eval(['oo_.bvar.forecast.with_shocks.Mean.' name ' = mean(sims, 2);']); diff --git a/matlab/bvar_irf.m b/matlab/bvar_irf.m index d742915d7..d7528f79b 100644 --- a/matlab/bvar_irf.m +++ b/matlab/bvar_irf.m @@ -138,9 +138,9 @@ save([ DirectoryName '/simulations.mat'], 'sampled_irfs'); % Save results in oo_ for i=1:ny - shock_name = options_.varobs(i, :); + shock_name = options_.varobs{i}; for j=1:ny - variable_name = options_.varobs(j, :); + variable_name = options_.varobs{j}; eval(['oo_.bvar.irf.Mean.' variable_name '.' shock_name ' = posterior_mean_irfs(' int2str(j) ',' int2str(i) ',:);']) eval(['oo_.bvar.irf.Median.' variable_name '.' shock_name ' = posterior_median_irfs(' int2str(j) ',' int2str(i) ',:);']) eval(['oo_.bvar.irf.Var.' variable_name '.' shock_name ' = posterior_variance_irfs(' int2str(j) ',' int2str(i) ',:);']) diff --git a/matlab/check_dsge_var_model.m b/matlab/check_dsge_var_model.m index b10815ac8..8b16ec09d 100644 --- a/matlab/check_dsge_var_model.m +++ b/matlab/check_dsge_var_model.m @@ -20,25 +20,25 @@ function check_dsge_var_model(Model, EstimatedParameters, BayesInfo) % along with Dynare. If not, see . if EstimatedParameters.nvn - error('DsgeVarLikelihood:: Measurement errors are not allowed!') + error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!') end if EstimatedParameters.ncn - error('DsgeVarLikelihood:: Measurement errors are not allowed!') + error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!') end if any(vec(Model.H)) - error('DsgeVarLikelihood:: Measurement errors are not allowed!') + error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!') end if EstimatedParameters.ncx - error('DsgeVarLikelihood:: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.') + error('Estimation::DsgeVarLikelihood: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.') end if Model.exo_nbr>1 && any(vec(tril(Model.Sigma_e,-1))) - error('DsgeVarLikelihood:: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.') + error('Estimation::DsgeVarLikelihood: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.') end if isequal(BayesInfo.with_trend,1) - error('DsgeVarLikelihood :: Linear trend is not yet implemented!') + error('Estimation::DsgeVarLikelihood: Linear trend is not yet implemented!') end diff --git a/matlab/check_file_extension.m b/matlab/check_file_extension.m index 63d627f0b..10d4c1c9b 100644 --- a/matlab/check_file_extension.m +++ b/matlab/check_file_extension.m @@ -19,6 +19,10 @@ function b = check_file_extension(file,type) % AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr +% Clean-up path +file = strrep(file, '../', ''); +file = strrep(file, './', ''); + remain = file; while ~isempty(remain) [ext, remain] = strtok(remain,'.'); diff --git a/matlab/check_list_of_variables.m b/matlab/check_list_of_variables.m index 48a3064a1..30805becb 100644 --- a/matlab/check_list_of_variables.m +++ b/matlab/check_list_of_variables.m @@ -46,7 +46,7 @@ if options_.dsge_var && options_.bayesian_irf msg = 1; end end - if size(varlist,1)~=size(options_.varobs) + if size(varlist,1)~=length(options_.varobs) msg = 1; end if msg @@ -55,7 +55,7 @@ if options_.dsge_var && options_.bayesian_irf skipline() end end - varlist = options_.varobs; + varlist = char(options_.varobs); return end @@ -128,7 +128,7 @@ elseif isempty(varlist) && isempty(options_.endo_vars_for_moment_computations_in if choice==1 varlist = M_.endo_names(1:M_.orig_endo_nbr, :); elseif choice==2 - varlist = options_.varobs; + varlist = char(options_.varobs); elseif choice==3 varlist = NaN; else @@ -138,15 +138,14 @@ elseif isempty(varlist) && isempty(options_.endo_vars_for_moment_computations_in end end end + if isnan(varlist) + edit([M_.fname '.mod']) + end + skipline() end - if isnan(varlist) - edit([M_.fname '.mod']) - end - skipline() end - function format_text(remain, max_number_of_words_per_line) index = 0; line_of_text = []; diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m index 5093ffc52..6482afb74 100644 --- a/matlab/compute_moments_varendo.m +++ b/matlab/compute_moments_varendo.m @@ -36,14 +36,14 @@ function oo_ = compute_moments_varendo(type,options_,M_,oo_,var_list_) if strcmpi(type,'posterior') posterior = 1; if nargin==4 - var_list_ = options_.varobs; + var_list_ = char(options_.varobs); end elseif strcmpi(type,'prior') posterior = 0; if nargin==4 var_list_ = options_.prior_analysis_endo_var_list; if isempty(var_list_) - options_.prior_analysis_var_list = options_.varobs; + options_.prior_analysis_var_list = char(options_.varobs); end end else diff --git a/matlab/display_estimation_results_table.m b/matlab/display_estimation_results_table.m index 0fecc8f38..334182bcb 100644 --- a/matlab/display_estimation_results_table.m +++ b/matlab/display_estimation_results_table.m @@ -105,7 +105,7 @@ if nvx disp(tit1) ip = nvx+1; for i=1:nvn - name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:)); + name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)}; if strcmp(field_name,'posterior') fprintf('%-*s %7.3f %8.4f %7.4f %4s %6.4f \n', ... header_width,name,bayestopt_.p1(ip), ... @@ -276,7 +276,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n'); ip = nvx+1; for i=1:nvn - idx = strmatch(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:),M_.endo_names); + idx = strmatch(options_.varobs{estim_params_.nvn_observable_correspondence(i,1)},M_.endo_names); fprintf(fidTeX,'$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',... deblank(M_.endo_names_tex(idx,:)), ... deblank(pnames(bayestopt_.pshape(ip)+1,:)), ... @@ -469,7 +469,7 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n'); ip = nvx+1; for i=1:nvn - idx = strmatch(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:),M_.endo_names); + idx = strmatch(options_.varobs{estim_params_.nvn_observable_correspondence(i,1)},M_.endo_names); fprintf(fidTeX,'$%s$ & %8.4f & %7.4f & %7.4f \\\\ \n',... deblank(M_.endo_names_tex(idx,:)), ... xparam1(ip),... diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m index e3ffc7cea..ee0c606a3 100644 --- a/matlab/dsge_likelihood.m +++ b/matlab/dsge_likelihood.m @@ -1,4 +1,4 @@ -function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info) +function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info) % Evaluates the posterior kernel of a dsge model using the specified % kalman_algo; the resulting posterior includes the 2*pi constant of the % likelihood function @@ -295,7 +295,7 @@ BayesInfo.mf = BayesInfo.mf1; % Define the constant vector of the measurement equation. if DynareOptions.noconstant - constant = zeros(DynareDataset.info.nvobs,1); + constant = zeros(DynareDataset.vobs,1); else if DynareOptions.loglinear constant = log(SteadyState(BayesInfo.mfys)); @@ -306,28 +306,28 @@ end % Define the deterministic linear trend of the measurement equation. if BayesInfo.with_trend - trend_coeff = zeros(DynareDataset.info.nvobs,1); + trend_coeff = zeros(DynareDataset.vobs,1); t = DynareOptions.trend_coeffs; for i=1:length(t) if ~isempty(t{i}) trend_coeff(i) = evalin('base',t{i}); end end - trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs]; + trend = repmat(constant,1,DynareDataset.nobs)+trend_coeff*[1:DynareDataset.nobs]; else - trend = repmat(constant,1,DynareDataset.info.ntobs); + trend = repmat(constant,1,DynareDataset.nobs); end % Get needed informations for kalman filter routines. start = DynareOptions.presample+1; -Z = BayesInfo.mf; % old mf -no_missing_data_flag = ~DynareDataset.missing.state; -mm = length(T); % old np -pp = DynareDataset.info.nvobs; +Z = BayesInfo.mf; +no_missing_data_flag = ~DatasetInfo.missing.state; +mm = length(T); +pp = DynareDataset.vobs; rr = length(Q); kalman_tol = DynareOptions.kalman_tol; riccati_tol = DynareOptions.riccati_tol; -Y = DynareDataset.data-trend; +Y = transpose(DynareDataset.data)-trend; %------------------------------------------------------------------------------ % 3. Initial condition of the Kalman filter @@ -406,7 +406,7 @@ switch DynareOptions.lik_init kalman_tol, riccati_tol, DynareOptions.presample, ... T,R,Q,H,Z,mm,pp,rr); else - [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ... + [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ... Y, 1, size(Y,2), ... zeros(mm,1), Pinf, Pstar, ... kalman_tol, riccati_tol, DynareOptions.presample, ... @@ -441,9 +441,9 @@ switch DynareOptions.lik_init % no need to test again for correlation elements correlated_errors_have_been_checked = 1; - [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,... - DynareDataset.missing.number_of_observations,... - DynareDataset.missing.no_more_missing_observations, ... + [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DatasetInfo.missing.aindex,... + DatasetInfo.missing.number_of_observations,... + DatasetInfo.missing.no_more_missing_observations, ... Y, 1, size(Y,2), ... zeros(mmm,1), Pinf, Pstar, ... kalman_tol, riccati_tol, DynareOptions.presample, ... @@ -523,7 +523,7 @@ if analytic_derivation, DLIK = []; AHess = []; iv = DynareResults.dr.restrict_var_list; - if nargin<8 || isempty(derivatives_info) + if nargin<9 || isempty(derivatives_info) [A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults); if ~isempty(EstimatedParameters.var_exo) indexo=EstimatedParameters.var_exo(:,1); @@ -655,10 +655,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter end else if 0 %DynareOptions.block - [err, LIK,lik] = block_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,... + [err, LIK,lik] = block_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,... T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar); else - [LIK,lik] = missing_observations_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... + [LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... a, Pstar, ... kalman_tol, DynareOptions.riccati_tol, ... DynareOptions.presample, ... @@ -733,7 +733,7 @@ if (kalman_algo==2) || (kalman_algo==4) end end - [LIK, lik] = univariate_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... + [LIK, lik] = univariate_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... a,Pstar, ... DynareOptions.kalman_tol, ... DynareOptions.riccati_tol, ... diff --git a/matlab/dsge_var_likelihood.m b/matlab/dsge_var_likelihood.m index 6c37fd1a4..f4cf0a8b2 100644 --- a/matlab/dsge_var_likelihood.m +++ b/matlab/dsge_var_likelihood.m @@ -1,4 +1,4 @@ -function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) +function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) % Evaluates the posterior kernel of the bvar-dsge model. % % INPUTS @@ -55,7 +55,11 @@ end nx = EstimatedParameters.nvx + EstimatedParameters.np; % Get the number of observed variables in the VAR model. -NumberOfObservedVariables = DynareDataset.info.nvobs; +NumberOfObservedVariables = DynareDataset.vobs; + +% Get the number of observations. +NumberOfObservations = DynareDataset.nobs; + % Get the number of lags in the VAR model. NumberOfLags = DynareOptions.dsge_varlag; @@ -110,8 +114,8 @@ Model.Sigma_e = Q; dsge_prior_weight = Model.params(dsge_prior_weight_idx); % Is the dsge prior proper? -if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.info.ntobs; - fval = objective_function_penalty_base+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables)); +if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/NumberOfObservations; + fval = objective_function_penalty_base+abs(NumberOfObservations*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables)); exit_flag = 0; info = 51; return @@ -197,9 +201,9 @@ assignin('base','GXX',GXX); assignin('base','GYX',GYX); if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is finite. - tmp0 = dsge_prior_weight*DynareDataset.info.ntobs*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ; - tmp1 = dsge_prior_weight*DynareDataset.info.ntobs*GYX + mYX; - tmp2 = inv(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX); + tmp0 = dsge_prior_weight*NumberOfObservations*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ; + tmp1 = dsge_prior_weight*NumberOfObservations*GYX + mYX; + tmp2 = inv(dsge_prior_weight*NumberOfObservations*GXX+mXX); SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0'); [SIGMAu_is_positive_definite, penalty] = ispd(SIGMAu); if ~SIGMAu_is_positive_definite @@ -208,28 +212,28 @@ if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model exit_flag = 0; return; end - SIGMAu = SIGMAu / (DynareDataset.info.ntobs*(1+dsge_prior_weight)); + SIGMAu = SIGMAu / (NumberOfObservations*(1+dsge_prior_weight)); PHI = tmp2*tmp1'; clear('tmp1'); - prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*DynareDataset.info.ntobs- ... + prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*NumberOfObservations- ... NumberOfObservedVariables*NumberOfLags ... +1-(1:NumberOfObservedVariables)'))); - prodlng2 = sum(gammaln(.5*(dsge_prior_weight*DynareDataset.info.ntobs- ... + prodlng2 = sum(gammaln(.5*(dsge_prior_weight*NumberOfObservations- ... NumberOfObservedVariables*NumberOfLags ... +1-(1:NumberOfObservedVariables)'))); - lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX)) ... - + .5*((dsge_prior_weight+1)*DynareDataset.info.ntobs-NumberOfParameters)*log(det((dsge_prior_weight+1)*DynareDataset.info.ntobs*SIGMAu)) ... - - .5*NumberOfObservedVariables*log(det(dsge_prior_weight*DynareDataset.info.ntobs*GXX)) ... - - .5*(dsge_prior_weight*DynareDataset.info.ntobs-NumberOfParameters)*log(det(dsge_prior_weight*DynareDataset.info.ntobs*(GYY-GYX*inv(GXX)*GYX'))) ... - + .5*NumberOfObservedVariables*DynareDataset.info.ntobs*log(2*pi) ... - - .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*DynareDataset.info.ntobs-NumberOfParameters) ... - + .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*DynareDataset.info.ntobs-NumberOfParameters) ... + lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX+mXX)) ... + + .5*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters)*log(det((dsge_prior_weight+1)*NumberOfObservations*SIGMAu)) ... + - .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX)) ... + - .5*(dsge_prior_weight*NumberOfObservations-NumberOfParameters)*log(det(dsge_prior_weight*NumberOfObservations*(GYY-GYX*inv(GXX)*GYX'))) ... + + .5*NumberOfObservedVariables*NumberOfObservations*log(2*pi) ... + - .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters) ... + + .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*NumberOfObservations-NumberOfParameters) ... - prodlng1 + prodlng2; else% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is infinite. iGXX = inv(GXX); SIGMAu = GYY - GYX*iGXX*transpose(GYX); PHI = iGXX*transpose(GYX); - lik = DynareDataset.info.ntobs * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ... - trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/DynareDataset.info.ntobs)); + lik = NumberOfObservations * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ... + trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/NumberOfObservations)); lik = .5*lik;% Minus likelihood end @@ -254,7 +258,7 @@ if (nargout==9) iGXX = inv(GXX); prior.SIGMAstar = GYY - GYX*iGXX*GYX'; prior.PHIstar = iGXX*transpose(GYX); - prior.ArtificialSampleSize = fix(dsge_prior_weight*DynareDataset.info.ntobs); + prior.ArtificialSampleSize = fix(dsge_prior_weight*NumberOfObservations); prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables; prior.iGXX = iGXX; end \ No newline at end of file diff --git a/matlab/dsgevar_posterior_density.m b/matlab/dsgevar_posterior_density.m index eca70f743..8e359f3ce 100644 --- a/matlab/dsgevar_posterior_density.m +++ b/matlab/dsgevar_posterior_density.m @@ -1,4 +1,4 @@ -function bvar = dsgevar_posterior_density(deep,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) +function bvar = dsgevar_posterior_density(deep,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) % This function characterizes the posterior distribution of a bvar with % a dsge prior (as in Del Negro and Schorfheide 2003) for a given value % of the deep parameters (structural parameters + the size of the @@ -38,7 +38,7 @@ dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight)); bvar.NumberOfLags = options_.varlag; -bvar.NumberOfVariables = size(options_.varobs,1); +bvar.NumberOfVariables = length(options_.varobs); bvar.Constant = 'no'; bvar.NumberOfEstimatedParameters = bvar.NumberOfLags*bvar.NumberOfVariables; if ~options_.noconstant @@ -47,7 +47,7 @@ if ~options_.noconstant bvar.NumberOfVariables; end -[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(deep',DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); +[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(deep',DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); % Conditionnal posterior density of the lagged matrices (given Sigma) -> % Matric-variate normal distribution. diff --git a/matlab/dyn_forecast.m b/matlab/dyn_forecast.m index bea68db53..5df0edbf8 100644 --- a/matlab/dyn_forecast.m +++ b/matlab/dyn_forecast.m @@ -81,8 +81,8 @@ switch task order_var = oo_.dr.order_var; i_var_obs = []; trend_coeffs = []; - for i=1:size(var_obs,1) - tmp = strmatch(var_obs(i,:),endo_names(i_var,:),'exact'); + for i=1:length(var_obs) + tmp = strmatch(var_obs{i},endo_names(i_var,:),'exact'); if ~isempty(tmp) i_var_obs = [ i_var_obs; tmp]; trend_coeffs = [trend_coeffs; oo_.Smoother.TrendCoeffs(i)]; diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m index 52d42aa0f..ed2e3d84a 100644 --- a/matlab/dynare_estimation.m +++ b/matlab/dynare_estimation.m @@ -105,7 +105,7 @@ if nnobs > 1 && horizon > 0 end endo_names = M_.endo_names; - n_varobs = size(options_.varobs,1); + n_varobs = length(options_.varobs); if isempty(var_list) var_list = endo_names; @@ -125,7 +125,7 @@ if nnobs > 1 && horizon > 0 IdObs = zeros(n_varobs,1); for j=1:n_varobs - iobs = strmatch(options_.varobs(j,:),var_list,'exact'); + iobs = strmatch(options_.varobs{j},var_list,'exact'); if ~isempty(iobs) IdObs(j,1) = iobs; end diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index df152bc30..a84f4a7fd 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -29,7 +29,7 @@ function dynare_estimation_1(var_list_,dname) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global M_ options_ oo_ estim_params_ bayestopt_ dataset_ +global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info % Set particle filter flag. if options_.order > 1 @@ -78,10 +78,20 @@ else objective_function = str2func('dsge_var_likelihood'); end -[dataset_,xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_); +[dataset_,dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_] = ... + dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_); if options_.dsge_var check_dsge_var_model(M_, estim_params_, bayestopt_); + if dataset_info.missing.state + error('Estimation::DsgeVarLikelihood: I cannot estimate a DSGE-VAR model with missing observations!') + end + if options_.noconstant + var_sample_moments(options_.dsge_varlag, -1, dataset_); + else + % The steady state is non zero ==> a constant in the VAR is needed! + var_sample_moments(options_.dsge_varlag, 0, dataset_); + end end %check for calibrated covariances before updating parameters @@ -114,14 +124,14 @@ if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(esti end data = dataset_.data; -rawdata = dataset_.rawdata; -data_index = dataset_.missing.aindex; -missing_value = dataset_.missing.state; +rawdata = dataset_.data; +data_index = dataset_info.missing.aindex; +missing_value = dataset_info.missing.state; % Set number of observations -gend = options_.nobs; +gend = dataset_.nobs; % Set the number of observed variables. -n_varobs = size(options_.varobs,1); +n_varobs = length(options_.varobs); % Get the number of parameters to be estimated. nvx = estim_params_.nvx; % Variance of the structural innovations (number of parameters). nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters). @@ -154,29 +164,10 @@ if ~isempty(estim_params_) M_ = set_all_parameters(xparam1,estim_params_,M_); end -% compute sample moments if needed (bvar-dsge) -if options_.dsge_var - if dataset_.missing.state - error('I cannot estimate a DSGE-VAR model with missing observations!') - end - if options_.noconstant - evalin('base',... - ['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ... - 'var_sample_moments(options_.first_obs,' ... - 'options_.first_obs+options_.nobs-1,options_.dsge_varlag,-1,' ... - 'options_.datafile, options_.varobs,options_.xls_sheet,options_.xls_range);']) - else% The steady state is non zero ==> a constant in the VAR is needed! - evalin('base',['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ... - 'var_sample_moments(options_.first_obs,' ... - 'options_.first_obs+options_.nobs-1,options_.dsge_varlag,0,' ... - 'options_.datafile, options_.varobs,options_.xls_sheet,options_.xls_range);']) - end -end - %% perform initial estimation checks; try - oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,M_,estim_params_,options_,bayestopt_,oo_); + oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,oo_); catch % if check fails, provide info on using calibration if present e = lasterror(); if full_calibration_detected %calibrated model present and no explicit starting values @@ -246,7 +237,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7); end [xparam1,fval,exitflag,output,lamdba,grad,hessian_fmincon] = ... - fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); case 2 error('ESTIMATION: mode_compute=2 option (Lester Ingber''s Adaptive Simulated Annealing) is no longer available') case 3 @@ -264,10 +255,10 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation optim_options = optimset(optim_options,'GradObj','on'); end if ~isoctave - [xparam1,fval,exitflag] = fminunc(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1,fval,exitflag] = fminunc(objective_function,xparam1,optim_options,dataset_,dataset_info_,options_,M_,estim_params_,bayestopt_,oo_); else % Under Octave, use a wrapper, since fminunc() does not have a 4th arg - func = @(x) objective_function(x, dataset_,options_,M_,estim_params_,bayestopt_,oo_); + func = @(x) objective_function(x, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); [xparam1,fval,exitflag] = fminunc(func,xparam1,optim_options); end case 4 @@ -306,7 +297,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation end % Call csminwell. [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ... - csminwel1(objective_function, xparam1, H0, analytic_grad, crit, nit, numgrad, epsilon, dataset_, options_, M_, estim_params_, bayestopt_, oo_); + csminwel1(objective_function, xparam1, H0, analytic_grad, crit, nit, numgrad, epsilon, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_); % Disp value at the mode. disp(sprintf('Objective function at mode: %f',fval)) case 5 @@ -334,7 +325,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation else nit=1000; end - [xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,oo_); if options_.analytic_derivation, options_.analytic_derivation = ana_deriv; end @@ -384,7 +375,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation end end % Evaluate the objective function. - fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + fval = feval(objective_function,xparam1,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,oo_); OldMode = fval; if ~exist('MeanPar','var') MeanPar = xparam1; @@ -416,8 +407,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation flag = 'LastCall'; end [xparam1,PostVar,Scale,PostMean] = ... - gmhmaxlik(objective_function,xparam1,[lb ub],gmhmaxlikOptions,Scale,flag,MeanPar,CovJump,dataset_,options_,M_,estim_params_,bayestopt_,oo_); - fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + gmhmaxlik(objective_function,xparam1,[lb ub],gmhmaxlikOptions,Scale,flag,MeanPar,CovJump,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); + fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); options_.mh_jscale = Scale; mouvement = max(max(abs(PostVar-OldPostVar))); skipline() @@ -436,11 +427,11 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation end [xparam1,PostVar,Scale,PostMean] = ... gmhmaxlik(objective_function,xparam1,[lb ub],... - gmhmaxlikOptions,Scale,flag,PostMean,PostVar,dataset_,options_,M_,estim_params_,bayestopt_,oo_); - fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + gmhmaxlikOptions,Scale,flag,PostMean,PostVar,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); + fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); options_.mh_jscale = Scale; mouvement = max(max(abs(PostVar-OldPostVar))); - fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); skipline() disp('========================================================== ') disp([' Change in the covariance matrix = ' num2str(mouvement) '.']) @@ -473,7 +464,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation if isfield(options_,'optim_opt') eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']); end - [xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); case 8 % Dynare implementation of the simplex algorithm. simplexOptions = options_.simplex; @@ -498,7 +489,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation end end end - [xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); case 9 % Set defaults H0 = 1e-4*ones(nx,1); @@ -523,7 +514,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation end warning('off','CMAES:NonfinitenessRange'); warning('off','CMAES:InitialSigma'); - [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,cmaesOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,cmaesOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); xparam1=BESTEVER.x; disp(sprintf('\n Objective function at mode: %f',fval)) case 10 @@ -554,16 +545,16 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation end simpsaOptionsList = options2cell(simpsaOptions); simpsaOptions = simpsaset(simpsaOptionsList{:}); - [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,simpsaOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,simpsaOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); case 11 options_.cova_compute = 0 ; - [xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_) ; + [xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) ; case 101 myoptions=soptions; myoptions(2)=1e-6; %accuracy of argument myoptions(3)=1e-6; %accuracy of function (see Solvopt p.29) myoptions(5)= 1.0; - [xparam1,fval]=solvopt(xparam1,objective_function,[],myoptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1,fval]=solvopt(xparam1,objective_function,[],myoptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); case 102 %simulating annealing % LB=zeros(size(xparam1))-20; @@ -600,10 +591,10 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation disp(['c vector ' num2str(c')]); [xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(objective_function,xparam1,maxy,rt_,epsilon,ns,nt ... - ,neps,maxevl,LB,UB,c,idisp ,t,vm,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + ,neps,maxevl,LB,UB,c,idisp ,t,vm,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); otherwise if ischar(options_.mode_compute) - [xparam1, fval, retcode ] = feval(options_.mode_compute,objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1, fval, retcode ] = feval(options_.mode_compute,objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); else error(['dynare_estimation:: mode_compute = ' int2str(options_.mode_compute) ' option is unknown!']) end @@ -614,11 +605,11 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation ana_deriv = options_.analytic_derivation; options_.analytic_derivation = 2; [junk1, junk2, hh] = feval(objective_function,xparam1, ... - dataset_,options_,M_,estim_params_,bayestopt_,oo_); + dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); options_.analytic_derivation = ana_deriv; else hh = reshape(hessian(objective_function,xparam1, ... - options_.gstep,dataset_,options_,M_,estim_params_,bayestopt_,oo_),nx,nx); + options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_),nx,nx); end end end @@ -698,7 +689,7 @@ end if options_.mode_check.status == 1 && ~options_.mh_posterior_mode_estimation ana_deriv = options_.analytic_derivation; options_.analytic_derivation = 0; - mode_check(objective_function,xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + mode_check(objective_function,xparam1,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); options_.analytic_derivation = ana_deriv; end @@ -725,14 +716,14 @@ if ~options_.cova_compute end if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation - %% display results table and store parameter estimates and standard errors in results + % display results table and store parameter estimates and standard errors in results oo_=display_estimation_results_table(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_,pnames,'Posterior','posterior'); - %% Laplace approximation to the marginal log density: + % Laplace approximation to the marginal log density: if options_.cova_compute estim_params_nbr = size(xparam1,1); scale_factor = -sum(log10(diag(invhess))); log_det_invhess = -estim_params_nbr*log(scale_factor)+log(det(scale_factor*invhess)); - likelihood = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood; skipline() disp(sprintf('Log data density [Laplace approximation] is %f.',oo_.MarginalDensity.LaplaceApproximation)) @@ -779,7 +770,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... ana_deriv = options_.analytic_derivation; options_.analytic_derivation = 0; if options_.cova_compute - feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); else error('I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.') end @@ -827,7 +818,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... if error_flag error('Estimation::mcmc: I cannot compute the posterior statistics!') end - prior_posterior_statistics('posterior',dataset_); + prior_posterior_statistics('posterior',dataset_,dataset_info); end xparam = get_posterior_parameters('mean'); M_ = set_all_parameters(xparam,estim_params_,M_); @@ -842,7 +833,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha > 0) && options_.load_mh_file)) ... || ~options_.smoother ) && options_.partial_information == 0 % to be fixed %% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable - [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,dataset_.info.ntobs,dataset_.data,dataset_.missing.aindex,dataset_.missing.state); + [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state); oo_.Smoother.SteadyState = ys; oo_.Smoother.TrendCoeffs = trend_coeff; oo_.Smoother.Variance = P; @@ -949,8 +940,11 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha %% Smooth observational errors... %% if options_.prefilter == 1 - yf = atT(bayestopt_.mf,:)+repmat(dataset_.descriptive.mean',1,gend); + yf = atT(bayestopt_.mf,:)+repmat(dataset_info.descriptive.mean',1,gend); elseif options_.loglinear == 1 + gend + bayestopt_.mfys + ys yf = atT(bayestopt_.mf,:)+repmat(log(ys(bayestopt_.mfys)),1,gend)+... trend_coeff*[1:gend]; else @@ -965,8 +959,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha number_of_plots_to_draw = number_of_plots_to_draw + 1; index = cat(1,index,i); end - eval(['oo_.SmoothedMeasurementErrors.' deblank(options_.varobs(i,:)) ... - ' = measurement_error(i,:)'';']); + eval(['oo_.SmoothedMeasurementErrors.' options_.varobs{i} ' = measurement_error(i,:)'';']); end if ~options_.nograph [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw); @@ -995,7 +988,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha hold on plot(1:gend,measurement_error(index(k),:),marker_string{2,1},'linewidth',1) hold off - name = deblank(options_.varobs(index(k),:)); + name = options_.varobs{index(k)}; if gend>1 xlim([1 gend]) end @@ -1009,7 +1002,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha set(gca,'XTickLabel',options_.XTickLabel) end if options_.TeX - idx = strmatch(options_.varobs(index(k),:),M_.endo_names,'exact'); + idx = strmatch(options_.varobs{index(k)},M_.endo_names,'exact'); texname = M_.endo_names_tex(idx,:); if isempty(TeXNAMES) TeXNAMES = ['$ ' deblank(texname) ' $']; @@ -1070,7 +1063,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha hold on plot(1:gend,rawdata(:,k),marker_string{2,1},'linewidth',1) hold off - name = deblank(options_.varobs(k,:)); + name = options_.varobs{k}; if isempty(NAMES) NAMES = name; else @@ -1084,7 +1077,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha xlim([1 gend]) end if options_.TeX - idx = strmatch(options_.varobs(k,:),M_.endo_names,'exact'); + idx = strmatch(options_.varobs{k},M_.endo_names,'exact'); texname = M_.endo_names_tex(idx,:); if isempty(TeXNAMES) TeXNAMES = ['$ ' deblank(texname) ' $']; diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index 5359470ba..874f96422 100644 --- a/matlab/dynare_estimation_init.m +++ b/matlab/dynare_estimation_init.m @@ -1,4 +1,4 @@ -function [dataset_, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, fake] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_) +function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, fake] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_) % function dynare_estimation_init(var_list_, gsa_flag) % preforms initialization tasks before estimation or @@ -41,47 +41,61 @@ hh = []; if isempty(gsa_flag) gsa_flag = 0; -else% Decide if a DSGE or DSGE-VAR has to be estimated. +else + % Decide if a DSGE or DSGE-VAR has to be estimated. if ~isempty(strmatch('dsge_prior_weight',M_.param_names)) options_.dsge_var = 1; end + % Get the list of the endogenous variables for which posterior statistics wil be computed var_list_ = check_list_of_variables(options_, M_, var_list_); options_.varlist = var_list_; end -% Get the indices of the observed variables in M_.endo_names. -options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1); -for i = 1:size(M_.endo_names,1) - tmp = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact'); - if ~isempty(tmp) - if length(tmp)>1 - skipline() - error(['Multiple declarations of ' deblank(M_.endo_names(i,:)) ' as an observed variable is not allowed!']) - end - options_.lgyidx2varobs(i) = tmp; +% Set the number of observed variables. +options_.number_of_observed_variables = length(options_.varobs); + +% Test if observed variables are declared. +if ~options_.number_of_observed_variables + error('VAROBS is missing!') +end + +% Check that each declared observed variable is also an endogenous variable. +for i = 1:options_.number_of_observed_variables + id = strmatch(options_.varobs{i}, M_.endo_names, 'exact'); + if isempty(id) + error(['Unknown variable (' options_.varobs{i} ')!']) end end +% Check that a variable is not declared as observed more than once. +if ~isequal(options_.varobs,unique(options_.varobs)) + for i = 1:options_.number_of_observed_variables + if length(strmatch(options_.varobs{i},options_.varobs))>1 + error(['A variable cannot be declared as observed more than once (' options_.varobs{i} ')!']) + end + end +end + +% Check the perturbation order (nonlinear filters with third order perturbation, or higher order, are not yet implemented). if options_.order>2 error(['I cannot estimate a model with a ' int2str(options_.order) ' order approximation of the model!']) end -% Set options_.lik_init equal to 3 if diffuse filter is used or -% kalman_algo refers to a diffuse filter algorithm. -if (options_.diffuse_filter==1) || (options_.kalman_algo > 2) - if options_.lik_init == 2 - error(['options diffuse_filter, lik_init and/or kalman_algo have ' ... - 'contradictory settings']) +% Set options_.lik_init equal to 3 if diffuse filter is used or kalman_algo refers to a diffuse filter algorithm. +if isequal(options_.diffuse_filter,1) || (options_.kalman_algo>2) + if isequal(options_.lik_init,2) + error(['options diffuse_filter, lik_init and/or kalman_algo have contradictory settings']) else options_.lik_init = 3; end end % If options_.lik_init == 1 -% set by default options_.qz_criterium to 1-1e-6 -% and check options_.qz_criterium < 1-eps if options_.lik_init == 1 -% Else set by default options_.qz_criterium to 1+1e-6 -if options_.lik_init == 1 +% set by default options_.qz_criterium to 1-1e-6 +% and check options_.qz_criterium < 1-eps if options_.lik_init == 1 +% Else +% set by default options_.qz_criterium to 1+1e-6 +if isequal(options_.lik_init,1) if isempty(options_.qz_criterium) options_.qz_criterium = 1-1e-6; elseif options_.qz_criterium > 1-eps @@ -111,9 +125,6 @@ else M_.dname = dname; end -% Set the number of observed variables. -n_varobs = size(options_.varobs,1); - % Set priors over the estimated parameters. if ~isempty(estim_params_) [xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_); @@ -297,7 +308,10 @@ end % storing prior parameters in results oo_.prior.mean = bayestopt_.p1; +oo_.prior.mode = bayestopt_.p5; oo_.prior.variance = diag(bayestopt_.p2.^2); +oo_.prior.hyperparameters.first = bayestopt_.p6; +oo_.prior.hyperparameters.second = bayestopt_.p7; % Is there a linear trend in the measurement equation? if ~isfield(options_,'trend_coeffs') % No! @@ -307,7 +321,7 @@ else% Yes! bayestopt_.trend_coeff = {}; trend_coeffs = options_.trend_coeffs; nt = length(trend_coeffs); - for i=1:n_varobs + for i=1:options_.number_of_observed_variables if i > length(trend_coeffs) bayestopt_.trend_coeff{i} = '0'; else @@ -326,17 +340,12 @@ nstatic = M_.nstatic; % Number of static variables. npred = M_.nspred; % Number of predetermined variables. nspred = M_.nspred; % Number of predetermined variables in the state equation. -% Test if observed variables are declared. -if isempty(options_.varobs) - error('VAROBS is missing') -end - % Setting resticted state space (observed + predetermined variables) var_obs_index = []; k1 = []; -for i=1:n_varobs - var_obs_index = [var_obs_index; strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:),'exact')]; - k1 = [k1; strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')]; +for i=1:options_.number_of_observed_variables + var_obs_index = [var_obs_index; strmatch(options_.varobs{i},M_.endo_names(dr.order_var,:),'exact')]; + k1 = [k1; strmatch(options_.varobs{i},M_.endo_names, 'exact')]; end % Define union of observed and state variables @@ -359,10 +368,8 @@ k3 = []; k3p = []; if options_.selected_variables_only for i=1:size(var_list_,1) - k3 = [k3; strmatch(var_list_(i,:),M_.endo_names(dr.order_var,:), ... - 'exact')]; - k3p = [k3; strmatch(var_list_(i,:),M_.endo_names, ... - 'exact')]; + k3 = [k3; strmatch(var_list_(i,:),M_.endo_names(dr.order_var,:), 'exact')]; + k3p = [k3; strmatch(var_list_(i,:),M_.endo_names, 'exact')]; end else k3 = (1:M_.endo_nbr)'; @@ -383,14 +390,12 @@ if options_.block == 1 bayestopt_.mf2 = var_obs_index; bayestopt_.mfys = k1; oo_.dr.restrict_columns = [size(i_posA,1)+(1:size(M_.state_var,2))]; - [k2, i_posA, i_posB] = union(k3p, M_.state_var', 'rows'); bayestopt_.smoother_var_list = [k3p(i_posA); M_.state_var(sort(i_posB))']; [junk,junk,bayestopt_.smoother_saved_var_list] = intersect(k3p,bayestopt_.smoother_var_list(:)); [junk,ic] = intersect(bayestopt_.smoother_var_list,M_.state_var); bayestopt_.smoother_restrict_columns = ic; - [junk,bayestopt_.smoother_mf] = ismember(k1, ... - bayestopt_.smoother_var_list); + [junk,bayestopt_.smoother_mf] = ismember(k1, bayestopt_.smoother_var_list); else k2 = union(var_obs_index,[M_.nstatic+1:M_.nstatic+M_.nspred]', 'rows'); % Set restrict_state to postion of observed + state variables in expanded state vector. @@ -404,13 +409,11 @@ else bayestopt_.mfys = k1; [junk,ic] = intersect(k2,nstatic+(1:npred)'); oo_.dr.restrict_columns = [ic; length(k2)+(1:nspred-npred)']; - bayestopt_.smoother_var_list = union(k2,k3); [junk,junk,bayestopt_.smoother_saved_var_list] = intersect(k3,bayestopt_.smoother_var_list(:)); [junk,ic] = intersect(bayestopt_.smoother_var_list,nstatic+(1:npred)'); bayestopt_.smoother_restrict_columns = ic; - [junk,bayestopt_.smoother_mf] = ismember(var_obs_index, ... - bayestopt_.smoother_var_list); + [junk,bayestopt_.smoother_mf] = ismember(var_obs_index, bayestopt_.smoother_var_list); end; if options_.analytic_derivation, @@ -440,37 +443,12 @@ if options_.analytic_derivation, end end -% Test if the data file is declared. -if isempty(options_.datafile) - if gsa_flag - dataset_ = []; -% rawdata = []; -% data_info = []; - return - else - error('datafile option is missing') - end -end - % If jscale isn't specified for an estimated parameter, use global option options_.jscale, set to 0.2, by default. k = find(isnan(bayestopt_.jscale)); bayestopt_.jscale(k) = options_.mh_jscale; -% Load and transform data. -transformation = []; -if options_.loglinear && ~options_.logdata - transformation = @log; -end -xls.sheet = options_.xls_sheet; -xls.range = options_.xls_range; - -if ~isfield(options_,'nobs') - options_.nobs = []; -end - -dataset_ = initialize_dataset(options_.datafile,options_.varobs,options_.first_obs,options_.nobs,transformation,options_.prefilter,xls); - -options_.nobs = dataset_.info.ntobs; +% Build the dataset +[dataset_, dataset_info] = makedataset(options_, options_.dsge_var*options_.dsge_varlag, gsa_flag); % setting steadystate_check_flag option if options_.diffuse_filter @@ -479,6 +457,7 @@ else steadystate_check_flag = 1; end +% If the steady state of the observed variables is non zero, set noconstant equal 0 () M = M_; nvx = estim_params_.nvx; ncx = estim_params_.ncx; @@ -502,5 +481,4 @@ else disp('variables in the model block, or drop the prefiltering.') error('The option "prefilter" is inconsistent with the non-zero mean measurement equations in the model.') end -end - +end \ No newline at end of file diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index bec4abe11..69fde1307 100644 --- a/matlab/dynare_identification.m +++ b/matlab/dynare_identification.m @@ -128,11 +128,11 @@ options_ = set_default_option(options_,'datafile',''); options_.mode_compute = 0; options_.plot_priors = 0; options_.smoother=1; -[dataset_,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_); +[dataset_,dataset_info,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_); options_ident.analytic_derivation_mode = options_.analytic_derivation_mode; if isempty(dataset_), dataset_.info.ntobs = periods; - dataset_.info.nvobs = rows(options_.varobs); + dataset_.info.nvobs = length(options_.varobs); dataset_.info.varobs = options_.varobs; dataset_.rawdata = []; dataset_.missing.state = 0; @@ -145,17 +145,8 @@ if isempty(dataset_), dataset_.missing.no_more_missing_observations = 1; dataset_.descriptive.mean = []; dataset_.data = []; - -% data_info.gend = periods; -% data_info.data = []; -% data_info.data_index = []; -% data_info.number_of_observations = periods*size(options_.varobs,1); -% data_info.no_more_missing_observations = 0; -% data_info.missing_value = 0; end -% results = prior_sampler(0,M_,bayestopt_,options_,oo_); - if prior_exist if any(bayestopt_.pshape > 0) if options_ident.prior_range @@ -278,7 +269,7 @@ if iload <=0, disp('Testing current parameter values') end [idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ... - identification_analysis(params,indx,indexo,options_ident,dataset_, prior_exist, name_tex,1); + identification_analysis(params,indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1); if info(1)~=0, skipline() disp('----------- ') @@ -321,7 +312,7 @@ if iload <=0, kk=kk+1; params = prior_draw(); [idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ... - identification_analysis(params,indx,indexo,options_ident,dataset_, prior_exist, name_tex,1); + identification_analysis(params,indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex,1); end end if info(1) @@ -375,7 +366,7 @@ if iload <=0, params = prior_draw(); end [dum1, ideJ, ideH, ideGP, dum2 , info] = ... - identification_analysis(params,indx,indexo,options_MC,dataset_, prior_exist, name_tex,0); + identification_analysis(params,indx,indexo,options_MC,dataset_, dataset_info, prior_exist, name_tex,0); if iteration==0 && info(1)==0, MAX_tau = min(SampleSize,ceil(MaxNumberOfBytes/(size(ideH.siH,1)*nparam)/8)); stoH = zeros([size(ideH.siH,1),nparam,MAX_tau]); @@ -546,7 +537,7 @@ if SampleSize > 1, disp(['Testing ',tittxt, '. Press ENTER']), pause(5), if ~iload, [idehess_max, idemoments_max, idemodel_max, idelre_max, derivatives_info_max] = ... - identification_analysis(pdraws(jmax,:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1); + identification_analysis(pdraws(jmax,:),indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex,1); save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_max', 'idemoments_max','idemodel_max', 'idelre_max', 'jmax', '-append'); end disp_identification(pdraws(jmax,:), idemodel_max, idemoments_max, name,1); @@ -561,7 +552,7 @@ if SampleSize > 1, disp(['Testing ',tittxt, '. Press ENTER']), pause(5), if ~iload, [idehess_min, idemoments_min, idemodel_min, idelre_min, derivatives_info_min] = ... - identification_analysis(pdraws(jmin,:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1); + identification_analysis(pdraws(jmin,:),indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1); save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_min', 'idemoments_min','idemodel_min', 'idelre_min', 'jmin', '-append'); end disp_identification(pdraws(jmin,:), idemodel_min, idemoments_min, name,1); @@ -576,7 +567,7 @@ if SampleSize > 1, disp(['Testing ',tittxt, '. Press ENTER']), pause(5), if ~iload, [idehess_(j), idemoments_(j), idemodel_(j), idelre_(j), derivatives_info_(j)] = ... - identification_analysis(pdraws(jcrit(j),:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1); + identification_analysis(pdraws(jcrit(j),:),indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1); end disp_identification(pdraws(jcrit(j),:), idemodel_(j), idemoments_(j), name,1); close all, diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m index 623f80ec5..1d9eeef7f 100644 --- a/matlab/dynare_sensitivity.m +++ b/matlab/dynare_sensitivity.m @@ -84,7 +84,7 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse, options_.mode_compute = 0; options_.filtered_vars = 1; options_.plot_priors = 0; - [dataset_,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_); + [dataset_,dataset_info,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_); % computes a first linear solution to set up various variables else if isempty(options_.qz_criterium) @@ -154,7 +154,7 @@ options_gsa = set_default_option(options_gsa,'namlagendo',[]); options_gsa = set_default_option(options_gsa,'namexo',[]); % RMSE mapping options_gsa = set_default_option(options_gsa,'lik_only',0); -options_gsa = set_default_option(options_gsa,'var_rmse',options_.varobs); +options_gsa = set_default_option(options_gsa,'var_rmse',char(options_.varobs)); options_gsa = set_default_option(options_gsa,'pfilt_rmse',0.1); options_gsa = set_default_option(options_gsa,'istart_rmse',options_.presample+1); options_gsa = set_default_option(options_gsa,'alpha_rmse',0.001); @@ -343,7 +343,7 @@ if options_gsa.rmse, end end - prior_posterior_statistics('gsa',dataset_); + prior_posterior_statistics('gsa',dataset_, dataset_info); if options_.bayesian_irf PosteriorIRF('gsa'); end @@ -366,7 +366,7 @@ if options_gsa.rmse, end clear a; % filt_mc_(OutputDirectoryName,data_info); - filt_mc_(OutputDirectoryName,options_gsa,dataset_); + filt_mc_(OutputDirectoryName,options_gsa,dataset_,dataset_info); end options_.opt_gsa = options_gsa; @@ -405,9 +405,9 @@ if options_gsa.glue, Obs.data = data; Obs.time = [1:gend]; Obs.num = gend; - for j=1:size(options_.varobs,1) - Obs.name{j} = deblank(options_.varobs(j,:)); - vj=deblank(options_.varobs(j,:)); + for j=1:length(options_.varobs) + Obs.name{j} = options_.varobs{j}; + vj = options_.varobs{j}; jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact'); js = strmatch(vj,lgy_,'exact'); @@ -445,7 +445,7 @@ if options_gsa.glue, ismoo(j)=jxj; end - jsmoo = size(options_.varobs,1); + jsmoo = length(options_.varobs); for j=1:M_.endo_nbr, if ~ismember(j,ismoo), jsmoo=jsmoo+1; @@ -470,10 +470,10 @@ if options_gsa.glue, Exo(j).name = deblank(tit(j,:)); end if ~options_gsa.ppost - Lik(size(options_.varobs,1)+1).name = 'logpo'; - Lik(size(options_.varobs,1)+1).ini = 'yes'; - Lik(size(options_.varobs,1)+1).isam = 1; - Lik(size(options_.varobs,1)+1).data = -logpo2; + Lik(length(options_.varobs)+1).name = 'logpo'; + Lik(length(options_.varobs)+1).ini = 'yes'; + Lik(length(options_.varobs)+1).isam = 1; + Lik(length(options_.varobs)+1).data = -logpo2; end Sam.name = bayestopt_.name; Sam.dim = [size(x) 0]; diff --git a/matlab/evaluate_likelihood.m b/matlab/evaluate_likelihood.m index 76fe3b427..14abbd0d4 100644 --- a/matlab/evaluate_likelihood.m +++ b/matlab/evaluate_likelihood.m @@ -37,7 +37,7 @@ function [llik,parameters] = evaluate_likelihood(parameters) global options_ M_ bayestopt_ oo_ estim_params_ -persistent dataset +persistent dataset dataset_info if nargin==0 parameters = 'posterior mode'; @@ -67,22 +67,9 @@ if ischar(parameters) end if isempty(dataset) - % Load and transform data. - transformation = []; - if options_.loglinear && ~options_.logdata - transformation = @log; - end - xls.sheet = options_.xls_sheet; - xls.range = options_.xls_range; - - if ~isfield(options_,'nobs') - options_.nobs = []; - end - - dataset = initialize_dataset(options_.datafile,options_.varobs,options_.first_obs,options_.nobs,transformation,options_.prefilter,xls); + [dataset, dataset_info] = makedataset(options_); end -llik = -dsge_likelihood(parameters,dataset,options_,M_,estim_params_,bayestopt_,oo_); +llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); ldens = evaluate_prior(parameters); -llik = llik - ldens; - +llik = llik - ldens; \ No newline at end of file diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m index 6b640b7b4..76b8444c0 100644 --- a/matlab/evaluate_smoother.m +++ b/matlab/evaluate_smoother.m @@ -44,14 +44,14 @@ function oo_=evaluate_smoother(parameters,var_list) global options_ M_ bayestopt_ oo_ estim_params_ % estim_params_ may be emty -persistent dataset_ +persistent dataset_ dataset_info if ischar(parameters) && strcmp(parameters,'calibration') options_.smoother=1; end if isempty(dataset_) || isempty(bayestopt_) - [dataset_,xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list, M_.fname, [], M_, options_, oo_, estim_params_, bayestopt_); + [dataset_,dataset_info,xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list, M_.fname, [], M_, options_, oo_, estim_params_, bayestopt_); end if nargin==0 @@ -88,7 +88,7 @@ if ischar(parameters) end [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = ... - DsgeSmoother(parameters,dataset_.info.ntobs,dataset_.data,dataset_.missing.aindex,dataset_.missing.state); + DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state); oo_.Smoother.SteadyState = ys; oo_.Smoother.TrendCoeffs = trend_coeff; diff --git a/matlab/forecast_graphs.m b/matlab/forecast_graphs.m index f95b77541..739e229d6 100644 --- a/matlab/forecast_graphs.m +++ b/matlab/forecast_graphs.m @@ -24,10 +24,6 @@ nr = 3; exo_nbr = M_.exo_nbr; endo_names = M_.endo_names; fname = M_.fname; -% $$$ varobs = options_.varobs; -% $$$ y = oo_.SmoothedVariables; -% $$$ ys = oo_.dr.ys; -% $$$ gend = size(y,2); yf = oo_.forecast.Mean; hpdinf = oo_.forecast.HPDinf; hpdsup = oo_.forecast.HPDsup; @@ -44,13 +40,6 @@ for i = 1:size(var_list) end nvar = length(i_var); -% $$$ % build trend for smoothed variables if necessary -% $$$ trend = zeros(size(varobs,1),10); -% $$$ if isfield(oo_.Smoother,'TrendCoeffs') -% $$$ trend_coeffs = oo_.Smoother.TrendCoeffs; -% $$$ trend = trend_coeffs*(gend-9:gend); -% $$$ end - % create subdirectory /graphs if id doesn't exist if ~exist(fname, 'dir') mkdir('.',fname); diff --git a/matlab/get_file_extension.m b/matlab/get_file_extension.m new file mode 100644 index 000000000..56c140d83 --- /dev/null +++ b/matlab/get_file_extension.m @@ -0,0 +1,42 @@ +function ext = get_file_extension(file) + +% returns the extension of a file. +% +% INPUTS +% o file string, name of the file +% +% OUTPUTS +% o ext string, extension. +% +% REMARKS +% If the provided file name has no extension, the routine will return an empty array. + +% Copyright (C) 2013 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 . + +% Clean-up path +file = strrep(file, '../', ''); +file = strrep(file, './', ''); + +remain = file; +while ~isempty(remain) + [ext, remain] = strtok(remain,'.'); +end + +if strcmp(ext,file) + ext = []; +end \ No newline at end of file diff --git a/matlab/get_posterior_parameters.m b/matlab/get_posterior_parameters.m index 8c4a9c8eb..b8676c2d4 100644 --- a/matlab/get_posterior_parameters.m +++ b/matlab/get_posterior_parameters.m @@ -50,7 +50,7 @@ end for i=1:nvn k1 = estim_params_.nvn_observable_correspondence(i,1); - name1 = deblank(options_.varobs(k1,:)); + name1 = options_.varobs{k1}; xparam(m) = eval(['oo_.posterior_' type '.measurement_errors_std.' name1]); m = m+1; end @@ -69,8 +69,8 @@ end for i=1:ncn k1 = estim_params_.corrn_observable_correspondence(i,1); k2 = estim_params_.corrn_observable_correspondence(i,2); - name1 = deblank(options_.varobs(k1,:)); - name2 = deblank(options_.varobs(k2,:)); + name1 = options_.varobs{k1}; + name2 = options_.varobs{k2}; xparam(m) = eval(['oo_.posterior_' type '.measurement_errors_corr.' name1 '_' name2]); m = m+1; end diff --git a/matlab/get_the_name.m b/matlab/get_the_name.m index bd6172a22..9251dd2f9 100644 --- a/matlab/get_the_name.m +++ b/matlab/get_the_name.m @@ -73,7 +73,7 @@ if k <= nvx texnam = ['$ SE_{' tname '} $']; end elseif k <= (nvx+nvn) - vname = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(k-estim_params_.nvx,1),:)); + vname = options_.varobs{estim_params_.nvn_observable_correspondence(k-estim_params_.nvx,1)}; nam=['SE_EOBS_',vname]; if TeX tname = deblank(M_.endo_names_tex(estim_params_.var_endo(k-estim_params_.nvx,1),:)); diff --git a/matlab/get_variables_list.m b/matlab/get_variables_list.m index 853b1169f..5d6a93c84 100644 --- a/matlab/get_variables_list.m +++ b/matlab/get_variables_list.m @@ -35,7 +35,7 @@ function [ivar,vartan,options_] = get_variables_list(options_,M_) varlist = options_.varlist; if isempty(varlist) - varlist = options_.varobs; + varlist = char(options_.varobs); options_.varlist = varlist; end nvar = rows(varlist); diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index f9c9c10fa..e9ee6c88d 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -33,6 +33,7 @@ global oo_ M_ options_ estim_params_ bayestopt_ estimation_info ex0_ ys0_ estim_params_ = []; bayestopt_ = []; options_.datafile = ''; +options_.dataset = []; options_.verbosity = 1; options_.terminal_condition = 0; options_.rplottype = 0; @@ -349,11 +350,13 @@ estimation_info.structural_innovation_corr_prior_index = {}; estimation_info.structural_innovation_corr_options_index = {}; estimation_info.structural_innovation_corr.range_index = {}; options_.initial_period = dates(1,1); -options_.dataset.firstobs = options_.initial_period; -options_.dataset.lastobs = NaN; +options_.dataset.file = []; +options_.dataset.series = []; +options_.dataset.firstobs = dates(); +options_.dataset.lastobs = dates(); options_.dataset.nobs = NaN; -options_.dataset.xls_sheet = NaN; -options_.dataset.xls_range = NaN; +options_.dataset.xls_sheet = []; +options_.dataset.xls_range = []; options_.Harvey_scale_factor = 10; options_.MaxNumberOfBytes = 1e6; options_.MaximumNumberOfMegaBytes = 111; @@ -364,7 +367,8 @@ options_.bayesian_th_moments = 0; options_.diffuse_filter = 0; options_.filter_step_ahead = []; options_.filtered_vars = 0; -options_.first_obs = 1; +options_.first_obs = NaN; +options_.nobs = NaN; options_.kalman_algo = 0; options_.fast_kalman = 0; options_.kalman_tol = 1e-10; diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m index ed6729362..329c6919d 100644 --- a/matlab/gsa/filt_mc_.m +++ b/matlab/gsa/filt_mc_.m @@ -1,4 +1,4 @@ -function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_) +function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info) % function [rmse_MC, ixx] = filt_mc_(OutDir) % inputs (from opt_gsa structure) % vvarvecm = options_gsa_.var_rmse; @@ -120,10 +120,10 @@ if ~loadSA, ys_mean=steady_(M_,options_,oo_); end % eval(options_.datafile) - Y = dataset_.data; - gend = dataset_.info.ntobs; - data_index = dataset_.missing.aindex; - missing_value = dataset_.missing.state; + Y = transpose(dataset_.data); + gend = dataset_.nobs; + data_index = dataset_info.missing.aindex; + missing_value = dataset_info.missing.state; for jx=1:gend, data_indx(jx,data_index{jx})=true; end %stock_gend=data_info.gend; %stock_data = data_info.data; diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index c30ed8a19..ce39b30cc 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -1,4 +1,4 @@ -function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist, name_tex, init) +function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex, init) % function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist, name_tex, init) % given the parameter vector params, wraps all identification analyses % @@ -134,17 +134,18 @@ if info(1)==0, options_.noprint = 1; options_.order = 1; options_.SpectralDensity.trigger = 0; - options_.periods = data_info.info.ntobs+100; + options_.periods = dataset_.nobs+100; if options_.kalman_algo > 2, options_.kalman_algo = 1; end analytic_derivation = options_.analytic_derivation; options_.analytic_derivation = -2; - info = stoch_simul(options_.varobs); - data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end); + info = stoch_simul(char(options_.varobs)); + dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end),dataset_.dates(1),dataset_.names,dataset_.tex); + %data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end); % datax=data; derivatives_info.no_DLIK=1; - [fval,DLIK,AHess,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info); + [fval,DLIK,AHess,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info); % fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_); options_.analytic_derivation = analytic_derivation; AHess=-AHess; diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m index 2fa56a2ba..766394fbc 100644 --- a/matlab/imcforecast.m +++ b/matlab/imcforecast.m @@ -102,30 +102,13 @@ if estimated_model error('imcforecast:: The dimension of the vector of parameters doesn''t match the number of estimated parameters!') end end - set_parameters(xparam); - - % Load and transform data. - transformation = []; - if options_.loglinear && ~options_.logdata - transformation = @log; - end - xls.sheet = options_.xls_sheet; - xls.range = options_.xls_range; - - if ~isfield(options_,'nobs') - options_.nobs = []; - end - - dataset_ = initialize_dataset(options_.datafile,options_.varobs,options_.first_obs,options_.nobs,transformation,options_.prefilter,xls); - - data = dataset_.data; - data_index = dataset_.missing.aindex; - gend = options_.nobs; - missing_value = dataset_.missing.state; - + [dataset_,dataset_info] = makedataset(options_); + data = transpose(dataset_.data); + data_index = dataset_info.missing.aindex; + gend = dataset_.nobs; + missing_value = dataset_info.missing.state; [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam,gend,data,data_index,missing_value); - trend = repmat(ys,1,options_cond_fcst.periods+1); for i=1:M_.endo_nbr j = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact'); @@ -134,7 +117,6 @@ if estimated_model end end trend = trend(oo_.dr.order_var,:); - InitState(:,1) = atT(:,end); else InitState(:,1) = zeros(M_.endo_nbr,1); diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m index 9dee5461e..c4cbf5ab9 100644 --- a/matlab/initial_estimation_checks.m +++ b/matlab/initial_estimation_checks.m @@ -1,11 +1,12 @@ -function DynareResults = initial_estimation_checks(objective_function,xparam1,DynareDataset,Model,EstimatedParameters,DynareOptions,BayesInfo,DynareResults) +function DynareResults = initial_estimation_checks(objective_function,xparam1,DynareDataset,DatasetInfo,Model,EstimatedParameters,DynareOptions,BayesInfo,DynareResults) % function initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations) % Checks data (complex values, ML evaluation, initial values, BK conditions,..) % % INPUTS % objective_function [function handle] of the objective function % xparam1: [vector] of parameters to be estimated -% DynareDataset: [structure] storing the dataset information +% DynareDataset: [dseries] object storing the dataset +% DataSetInfo: [structure] storing informations about the sample. % Model: [structure] decribing the model % EstimatedParameters [structure] characterizing parameters to be estimated % DynareOptions [structure] describing the options @@ -35,11 +36,11 @@ function DynareResults = initial_estimation_checks(objective_function,xparam1,Dy % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -if DynareDataset.info.nvobs>Model.exo_nbr+EstimatedParameters.nvn +if DynareDataset.vobs>Model.exo_nbr+EstimatedParameters.nvn error(['initial_estimation_checks:: Estimation can''t take place because there are less declared shocks than observed variables!']) end - -if DynareDataset.info.nvobs>length(find(diag(Model.Sigma_e)))+EstimatedParameters.nvn + +if DynareDataset.vobs>length(find(diag(Model.Sigma_e)))+EstimatedParameters.nvn error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!']) end @@ -72,7 +73,7 @@ end % Evaluate the likelihood. ana_deriv = DynareOptions.analytic_derivation; DynareOptions.analytic_derivation=0; -[fval,junk1,junk2,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); +[fval,junk1,junk2,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); DynareOptions.analytic_derivation=ana_deriv; if DynareOptions.dsge_var || strcmp(func2str(objective_function),'non_linear_dsge_likelihood') diff --git a/matlab/load_m_file_data.m b/matlab/load_m_file_data.m index b69258583..20e1b28a7 100644 --- a/matlab/load_m_file_data.m +++ b/matlab/load_m_file_data.m @@ -36,7 +36,7 @@ function [freq,init,data,varlist,tex] = load_m_file_data(file) if isoctave run(file); else - [basename, ext] = strtok(file,'.'); + basename = file(1:end-2); run(basename); end diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m index 8c6ed5a9f..7fb056553 100644 --- a/matlab/metropolis_hastings_initialization.m +++ b/matlab/metropolis_hastings_initialization.m @@ -1,7 +1,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... - metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_) + metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) %function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = -% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_,options_,M_,estim_params_,bayestopt_,oo_) +% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_,dataset_info,,options_,M_,estim_params_,bayestopt_,oo_) % Metropolis-Hastings initialization. % % INPUTS @@ -110,7 +110,7 @@ if ~options_.load_mh_file && ~options_.mh_recover candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar); if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2)) ix2(j,:) = candidate; - ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,options_,M_,estim_params_,bayestopt_,oo_); + ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); if ~isfinite(ilogpo2(j)) % if returned log-density is % Inf or Nan (penalized value) validate = 0; @@ -152,7 +152,7 @@ if ~options_.load_mh_file && ~options_.mh_recover candidate = transpose(xparam1(:));% if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2)) ix2 = candidate; - ilogpo2 = - feval(TargetFun,ix2',dataset_,options_,M_,estim_params_,bayestopt_,oo_); + ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); disp('Estimation::mcmc: Initialization at the posterior mode.') skipline() fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']); diff --git a/matlab/mode_check.m b/matlab/mode_check.m index 786acb7da..2fa6f4a9d 100644 --- a/matlab/mode_check.m +++ b/matlab/mode_check.m @@ -1,4 +1,4 @@ -function mode_check(fun,x,hessian,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) +function mode_check(fun,x,hessian,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) % Checks the estimated ML mode or Posterior mode. %@info: @@ -62,7 +62,7 @@ if ~isempty(hessian); [ s_min, k ] = min(diag(hessian)); end -fval = feval(fun,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); +fval = feval(fun,x,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); if ~isempty(hessian); skipline() @@ -138,7 +138,7 @@ for plt = 1:nbplt, end for i=1:length(z) xx(kk) = z(i); - [fval, junk1, junk2, exit_flag] = feval(fun,xx,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); + [fval, junk1, junk2, exit_flag] = feval(fun,xx,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); if exit_flag y(i,1) = fval; else diff --git a/matlab/mr_gstep.m b/matlab/mr_gstep.m index 690ffa01b..82c8d70f4 100644 --- a/matlab/mr_gstep.m +++ b/matlab/mr_gstep.m @@ -1,9 +1,17 @@ -function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) +function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin) % function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin) % % Gibbs type step in optimisation +% +% varargin{1} --> DynareDataset +% varargin{2} --> DatasetInfo +% varargin{3} --> DynareOptions +% varargin{4} --> Model +% varargin{5} --> EstimatedParameters +% varargin{6} --> BayesInfo +% varargin{1} --> DynareResults -% Copyright (C) 2006-2012 Dynare Team +% Copyright (C) 2006-2014 Dynare Team % % This file is part of Dynare. % @@ -22,7 +30,7 @@ function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,DynareDataset,DynareOptions,Mod n=size(x,1); if isempty(h1), - h1=DynareOptions.gradient_epsilon*ones(n,1); + h1=varargin{3}.gradient_epsilon*ones(n,1); end @@ -31,7 +39,7 @@ if isempty(htol0) else htol = htol0; end -f0=feval(func0,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); +f0=feval(func0,x,varargin{:}); xh1=x; f1=zeros(size(f0,1),n); @@ -45,10 +53,10 @@ while i htol - [f0 x fc retcode] = csminit1(func0,x,f0,gg,0,diag(hh),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); + [f0 x fc retcode] = csminit1(func0,x,f0,gg,0,diag(hh),varargin{:}); ig(i)=1; - fprintf(['Done for param %s = %8.4f\n'],BayesInfo.name{i},x(i)) + fprintf(['Done for param %s = %8.4f\n'],varargin{6}.name{i},x(i)) end xh1=x; end diff --git a/matlab/mr_hessian.m b/matlab/mr_hessian.m index a428a3c92..1cdbb05b1 100644 --- a/matlab/mr_hessian.m +++ b/matlab/mr_hessian.m @@ -1,10 +1,10 @@ -function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) +function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,varargin) % [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,varargin) % % numerical gradient and Hessian, with 'automatic' check of numerical % error % -% adapted from Michel Juillard original rutine hessian.m +% adapted from Michel Juillard original routine hessian.m % % func = function handle. The function must give two outputs: % - the log-likelihood AND the single contributions at times t=1,...,T @@ -22,9 +22,17 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hf % htol0 = 'precision' of increment of function values for numerical % derivatives % -% varargin: other parameters of func +% varargin{1} --> DynareDataset +% varargin{2} --> DatasetInfo +% varargin{3} --> DynareOptions +% varargin{4} --> Model +% varargin{5} --> EstimatedParameters +% varargin{6} --> BayesInfo +% varargin{1} --> DynareResults -% Copyright (C) 2004-2012 Dynare Team + + +% Copyright (C) 2004-2014 Dynare Team % % This file is part of Dynare. % @@ -45,16 +53,16 @@ persistent h1 htol n=size(x,1); if init - gstep_=DynareOptions.gstep; + gstep_=varargin{3}.gstep; htol = 1.e-4; - h1=DynareOptions.gradient_epsilon*ones(n,1); + h1=varargin{3}.gradient_epsilon*ones(n,1); return end -[f0, ff0]=feval(func,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); -h2=BayesInfo.ub-BayesInfo.lb; -hmax=BayesInfo.ub-x; -hmax=min(hmax,x-BayesInfo.lb); +[f0, ff0]=feval(func,x,varargin{:}); +h2=varargin{6}.ub-varargin{6}.lb; +hmax=varargin{6}.ub-x; +hmax=min(hmax,x-varargin{6}.lb); if isempty(ff0), outer_product_gradient=0; else @@ -83,7 +91,7 @@ while i(3*htol)), + if (h1(i)<1.e-12*min(1,h2(i)) && h1(i)<0.5*hmax(i)) ic=1; hcheck=1; end @@ -141,7 +149,7 @@ while i. n_chains = length(options.ms.ms_chain); - nvars = size(options.varobs,1); + nvars = length(options.varobs); fh = fopen(fname,'w'); %/******************************************************************************/ diff --git a/matlab/ms-sbvar/msstart_setup.m b/matlab/ms-sbvar/msstart_setup.m index 2d322211d..37b6e27e9 100644 --- a/matlab/ms-sbvar/msstart_setup.m +++ b/matlab/ms-sbvar/msstart_setup.m @@ -48,16 +48,16 @@ end %1 CBO output gap -- log(x_t)-log(x_t potential) %2 GDP deflator -- (P_t/P_{t-1})^4-1.0 %2 FFR/100. -options_.ms.vlist = [1:size(options_.varobs,1)]; % 1: U; 4: PCE inflation. -options_.ms.varlist=cellstr(options_.varobs); -options_.ms.log_var = sort(varlist_indices(options_.ms.vlistlog,options_.varobs)); % subset of "options_.ms.vlist. Variables in log level so that differences are in **monthly** growth, unlike R and U which are in annual percent (divided by 100 already). +options_.ms.vlist = [1:length(options_.varobs)]; % 1: U; 4: PCE inflation. +options_.ms.varlist=cellstr(options_.varobs'); +options_.ms.log_var = sort(varlist_indices(options_.ms.vlistlog,char(options_.varobs))); % subset of "options_.ms.vlist. Variables in log level so that differences are in **monthly** growth, unlike R and U which are in annual percent (divided by 100 already). options_.ms.percent_var =setdiff(options_.ms.vlist,options_.ms.log_var); %options_.ms.restriction_fname='ftd_upperchol3v'; %Only used by msstart2.m. ylab = options_.ms.varlist; xlab = options_.ms.varlist; %---------------- -nvar = size(options_.varobs,1); % number of endogenous variables +nvar = length(options_.varobs); % number of endogenous variables nlogeno = length(options_.ms.log_var); % number of endogenous variables in options_.ms.log_var npereno = length(options_.ms.percent_var); % number of endogenous variables in options_.ms.percent_var if (nvar~=(nlogeno+npereno)) diff --git a/matlab/newrat.m b/matlab/newrat.m index 1df015776..424db0c1c 100644 --- a/matlab/newrat.m +++ b/matlab/newrat.m @@ -1,4 +1,4 @@ -function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ftol0, nit, flagg, DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) +function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ftol0, nit, flagg, varargin) % [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit, flagg, varargin) % % Optimiser with outer product gradient and with sequences of univariate steps @@ -21,9 +21,16 @@ function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ft % with correlation structure as from outer product gradient, % flagg = 2, full numerical Hessian % -% varargin = list of parameters for func0 +% varargin{1} --> DynareDataset +% varargin{2} --> DatasetInfo +% varargin{3} --> DynareOptions +% varargin{4} --> Model +% varargin{5} --> EstimatedParameters +% varargin{6} --> BayesInfo +% varargin{1} --> DynareResults -% Copyright (C) 2004-2013 Dynare Team + +% Copyright (C) 2004-2014 Dynare Team % % This file is part of Dynare. % @@ -52,19 +59,19 @@ ftol=ftol0; gtol=1.e-3; htol=htol_base; htol0=htol_base; -gibbstol=length(BayesInfo.pshape)/50; %25; +gibbstol=length(varargin{6}.pshape)/50; %25; % func0 = str2func([func2str(func0),'_hh']); % func0 = func0; -[fval0,gg,hh]=feval(func0,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); +[fval0,gg,hh]=feval(func0,x,varargin{:}); fval=fval0; % initialize mr_gstep and mr_hessian outer_product_gradient=1; if isempty(hh) - mr_hessian(1,x,[],[],[],DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); - [dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,flagit,htol,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); + mr_hessian(1,x,[],[],[],varargin{:}); + [dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,flagit,htol,varargin{:}); if isempty(dum), outer_product_gradient=0; igg = 1e-4*eye(nx); @@ -111,9 +118,9 @@ while norm(gg)>gtol && check==0 && jit1 disp('Gradient step!!') else @@ -122,28 +129,26 @@ while norm(gg)>gtol && check==0 && jit1 && (fval0(icount-1)-fval0(icount))>1) || ((fval0(icount)-fval)<1.e-2*(gg'*(H*gg))/2 && igibbs), - if length(find(ig))=ftol disp('Diagonal Hessian successful') end @@ -152,7 +157,7 @@ while norm(gg)>gtol && check==0 && jit=ftol disp('Gradient direction successful') end @@ -165,14 +170,14 @@ while norm(gg)>gtol && check==0 && jit0 - [dum, gg, htol0, igg, hhg,h1]=mr_hessian(0,xparam1,func0,flagg,ftol0,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); + [dum, gg, htol0, igg, hhg,h1]=mr_hessian(0,xparam1,func0,flagg,ftol0,varargin{:}); if flagg==2 hh = reshape(dum,nx,nx); ee=eig(hh); @@ -208,7 +213,7 @@ while norm(gg)>gtol && check==0 && jitgtol && check==0 && jit MAX_momentsno || b == B) - % stock = stock_moments(1:irun(8)-1); - % ifil(8) = ifil(8) + 1; - % save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock'); - % if RemoteFlag==1, - % OutputFileName_moments = [OutputFileName_moments; {[DirectoryName filesep], [M_.fname '_moments' int2str(ifil(8)) '.mat']}]; - % end - % irun(8) = 1; - % end - - % DirectoryName=TempPath; - - dyn_waitbar((b-fpar+1)/(B-fpar+1),h); - end myoutput.ifil=ifil; @@ -332,7 +313,6 @@ if RemoteFlag==1, OutputFileName_param; OutputFileName_forc_mean; OutputFileName_forc_point]; - % OutputFileName_moments]; end dyn_waitbar_close(h); diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m index f524b44c1..c2a49e7a2 100644 --- a/matlab/random_walk_metropolis_hastings.m +++ b/matlab/random_walk_metropolis_hastings.m @@ -1,4 +1,4 @@ -function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_) +function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) %function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_) % Random walk Metropolis-Hastings algorithm. % @@ -60,7 +60,7 @@ objective_function_penalty_base = Inf; % Initialization of the random walk metropolis-hastings chains. [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... - metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2); @@ -101,6 +101,7 @@ localVars = struct('TargetFun', TargetFun, ... 'InitSizeArray',InitSizeArray, ... 'record', record, ... 'dataset_', dataset_, ... + 'dataset_info', dataset_info, ... 'options_', options_, ... 'M_',M_, ... 'bayestopt_', bayestopt_, ... diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m index 481e92b77..62e1728eb 100644 --- a/matlab/random_walk_metropolis_hastings_core.m +++ b/matlab/random_walk_metropolis_hastings_core.m @@ -88,6 +88,7 @@ d=myinputs.d; InitSizeArray=myinputs.InitSizeArray; record=myinputs.record; dataset_ = myinputs.dataset_; +dataset_info = myinputs.dataset_info; bayestopt_ = myinputs.bayestopt_; estim_params_ = myinputs.estim_params_; options_ = myinputs.options_; @@ -164,7 +165,7 @@ for b = fblck:nblck, par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n); if all( par(:) > mh_bounds(:,1) ) && all( par(:) < mh_bounds(:,2) ) try - logpost = - feval(TargetFun, par(:),dataset_,options_,M_,estim_params_,bayestopt_,oo_); + logpost = - feval(TargetFun, par(:),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); catch logpost = -inf; end diff --git a/matlab/read_variables.m b/matlab/read_variables.m index 23d7fe635..6859a6795 100644 --- a/matlab/read_variables.m +++ b/matlab/read_variables.m @@ -42,7 +42,7 @@ if ~isempty(directory) end dyn_size_01 = size(dyn_data_01,1); -var_size_01 = size(var_names_01,1); +var_size_01 = length(var_names_01); % Auto-detect extension if not provided if isempty(extension) @@ -71,7 +71,7 @@ switch (extension) case '.m' eval(basename); for dyn_i_01=1:var_size_01 - dyn_tmp_01 = eval(var_names_01(dyn_i_01,:)); + dyn_tmp_01 = eval(var_names_01{dyn_i_01}); if length(dyn_tmp_01) > dyn_size_01 && dyn_size_01 > 0 cd(old_pwd) error('data size is too large') @@ -81,7 +81,7 @@ switch (extension) case '.mat' s = load(basename); for dyn_i_01=1:var_size_01 - dyn_tmp_01 = s.(deblank(var_names_01(dyn_i_01,:))); + dyn_tmp_01 = s.(var_names_01{dyn_i_01}); if length(dyn_tmp_01) > dyn_size_01 && dyn_size_01 > 0 cd(old_pwd) error('data size is too large') @@ -106,9 +106,8 @@ switch (extension) end case '.csv' [freq,init,data,varlist] = load_csv_file_data(fullname); - %var_names_01 = deblank(var_names_01); for dyn_i_01=1:var_size_01 - iv = strmatch(strtrim(var_names_01(dyn_i_01,:)),varlist,'exact'); + iv = strmatch(var_names_01{dyn_i_01},varlist,'exact'); if ~isempty(iv) dyn_tmp_01 = [data(:,iv)]'; if length(dyn_tmp_01) > dyn_size_01 && dyn_size_01 > 0 @@ -118,7 +117,7 @@ switch (extension) dyn_data_01(:,dyn_i_01) = dyn_tmp_01; else cd(old_pwd) - error([strtrim(var_names_01(dyn_i_01,:)) ' not found in ' fullname]) + error([var_names_01{dyn_i_01} ' not found in ' fullname]) end end otherwise diff --git a/matlab/set_default_option.m b/matlab/set_default_option.m index 3c8303f76..760d74fa8 100644 --- a/matlab/set_default_option.m +++ b/matlab/set_default_option.m @@ -33,6 +33,17 @@ function options=set_default_option(options,field,default) if ~isfield(options,field) options.(field) = default; + return +end + +if isempty(options.(field)) + options.(field) = default; + return +end + +if isnan(options.(field)) + options.(field) = default; + return end % 06/07/03 MJ added ; to eval expression \ No newline at end of file diff --git a/matlab/set_prior.m b/matlab/set_prior.m index c438b2d3b..253b51c24 100644 --- a/matlab/set_prior.m +++ b/matlab/set_prior.m @@ -76,12 +76,12 @@ end if nvn estim_params_.nvn_observable_correspondence=NaN(nvn,1); % stores number of corresponding observable if isequal(M_.H,0) %if no previously set measurement error, initialize H - nvarobs = size(options_.varobs,1); + nvarobs = length(options_.varobs); M_.H = zeros(nvarobs,nvarobs); M_.Correlation_matrix_ME = eye(nvarobs); end for i=1:nvn - obsi_ = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),deblank(options_.varobs),'exact'); + obsi_ = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),options_.varobs,'exact'); if isempty(obsi_) error(['The variable ' deblank(M_.endo_names(estim_params_.var_endo(i,1),:)) ' has to be declared as observable since you assume a measurement error on it.']) end @@ -96,7 +96,7 @@ if nvn bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_endo(:,8)]; %take generalized distribution into account bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_endo(:,9)]; %take generalized distribution into account bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_endo(:,10)]; - bayestopt_.name = [ bayestopt_.name; cellstr(options_.varobs(estim_params_.nvn_observable_correspondence,:))]; + bayestopt_.name = [ bayestopt_.name; transpose(options_.varobs(estim_params_.nvn_observable_correspondence))]; end if ncx xparam1 = [xparam1; estim_params_.corrx(:,3)]; @@ -115,7 +115,7 @@ end if ncn estim_params_.corrn_observable_correspondence=NaN(ncn,2); if isequal(M_.H,0) - nvarobs = size(options_.varobs,1); + nvarobs = length(options_.varobs); M_.H = zeros(nvarobs,nvarobs); M_.Correlation_matrix_ME = eye(nvarobs); end @@ -134,8 +134,8 @@ if ncn for i=1:ncn k1 = estim_params_.corrn(i,1); k2 = estim_params_.corrn(i,2); - obsi1 = strmatch(deblank(M_.endo_names(k1,:)),deblank(options_.varobs),'exact'); %find correspondence to varobs to construct H in set_all_paramters - obsi2 = strmatch(deblank(M_.endo_names(k2,:)),deblank(options_.varobs),'exact'); + obsi1 = strmatch(deblank(M_.endo_names(k1,:)),options_.varobs,'exact'); %find correspondence to varobs to construct H in set_all_paramters + obsi2 = strmatch(deblank(M_.endo_names(k2,:)),options_.varobs,'exact'); estim_params_.corrn_observable_correspondence(i,:)=[obsi1,obsi2]; %save correspondence end end diff --git a/matlab/simplex_optimization_routine.m b/matlab/simplex_optimization_routine.m index da073cb19..9d60b5774 100644 --- a/matlab/simplex_optimization_routine.m +++ b/matlab/simplex_optimization_routine.m @@ -11,8 +11,15 @@ function [x,fval,exitflag] = simplex_optimization_routine(objective_function,x,o % o objective_function [string] Name of the objective function to be minimized. % o x [double] n*1 vector, starting guess of the optimization routine. % o options [structure] Options of this implementation of the simplex algorithm. -% o varargin [cell of structures] Structures to be passed to the objective function: dataset_, -% options_, M_, estim_params_, bayestopt_, and oo_. +% o varargin [cell of structures] Structures to be passed to the objective function. +% +% varargin{1} --> DynareDataset +% varargin{2} --> DatasetInfo +% varargin{3} --> DynareOptions +% varargin{4} --> Model +% varargin{5} --> EstimatedParameters +% varargin{6} --> BayesInfo +% varargin{1} --> DynareResults % % OUTPUTS % o x [double] n*1 vector, estimate of the optimal inputs. @@ -187,7 +194,7 @@ else disp(['Current parameter values: ']) fprintf(1,'%s: \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \n','Names','Best point', 'Worst point', 'Mean values', 'Min values', 'Max values'); for i=1:number_of_variables - fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{5}.name{i},v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2)); + fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{6}.name{i},v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2)); end skipline() end @@ -399,7 +406,7 @@ while (func_count < max_func_calls) && (iter_count < max_iterations) && (simplex disp(['Current parameter values: ']) fprintf(1,'%s: \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \n','Names','Best point', 'Worst point', 'Mean values', 'Min values', 'Max values'); for i=1:number_of_variables - fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{5}.name{i}, v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2)); + fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{6}.name{i}, v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2)); end skipline() end @@ -425,7 +432,7 @@ while (func_count < max_func_calls) && (iter_count < max_iterations) && (simplex disp(['values for the control variables. ']) disp(['New value of delta (size of the new simplex) is: ']) for i=1:number_of_variables - fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{5}.name{i}, delta(i)); + fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{6}.name{i}, delta(i)); end end % Reset counters @@ -472,7 +479,7 @@ while (func_count < max_func_calls) && (iter_count < max_iterations) && (simplex disp(['values for the control variables. ']) disp(['New value of delta (size of the new simplex) is: ']) for i=1:number_of_variables - fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{5}.name{i}, delta(i)); + fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{6}.name{i}, delta(i)); end end % Reset counters diff --git a/matlab/simulated_moment_uncertainty.m b/matlab/simulated_moment_uncertainty.m index 0815c7a6b..d41f13835 100644 --- a/matlab/simulated_moment_uncertainty.m +++ b/matlab/simulated_moment_uncertainty.m @@ -30,7 +30,7 @@ for j=1:replic; options_.noprint = 1; options_.order = 1; options_.periods = periods; - info = stoch_simul(options_.varobs); + info = stoch_simul(char(options_.varobs)); dum=[oo_.mean; dyn_vech(oo_.var)]; sd = sqrt(diag(oo_.var)); for i=1:options_.ar; diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index 8b2144651..2009ca08d 100644 --- a/matlab/stoch_simul.m +++ b/matlab/stoch_simul.m @@ -115,7 +115,7 @@ if ~options_.noprint skipline() if isfield(options_,'varobs')&& ~isempty(options_.varobs) - PCL_varobs=options_.varobs; + PCL_varobs=char(options_.varobs); disp('OBSERVED VARIABLES') else PCL_varobs=M_.endo_names; diff --git a/matlab/subset.m b/matlab/subset.m index 8e419581c..dedbdbfba 100644 --- a/matlab/subset.m +++ b/matlab/subset.m @@ -20,8 +20,6 @@ function jndx = subset() global options_ estim_params_ M_ ExcludedParamNames = options_.ExcludedParams; -VarObs = options_.varobs; -VarExo = M_.exo_names; info = options_.ParamSubSet; nvx = estim_params_.nvx; diff --git a/matlab/utilities/dataset/compute_corr.m b/matlab/utilities/dataset/compute_corr.m deleted file mode 100644 index d6469baef..000000000 --- a/matlab/utilities/dataset/compute_corr.m +++ /dev/null @@ -1,59 +0,0 @@ -function dataset_ = compute_corr(dataset_) -% Computes the correlation matrix of the sample (possibly with missing observations). - -%@info: -%! @deftypefn {Function File} {@var{dataset_} =} compute_corr(@var{dataset_}) -%! @anchor{compute_corr} -%! This function computes covariance matrix of the sample (possibly with missing observations). -%! -%! @strong{Inputs} -%! @table @var -%! @item dataset_ -%! Dynare structure describing the dataset, built by @ref{initialize_dataset} -%! @end table -%! -%! @strong{Outputs} -%! @table @var -%! @item dataset_ -%! Dynare structure describing the dataset, built by @ref{initialize_dataset} -%! @end table -%! -%! @strong{This function is called by:} -%! @ref{descriptive_statistics}. -%! -%! @strong{This function calls:} -%! @ref{ndim}, @ref{compute_cova}. -%! -%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.corr} is a -%! @tex{n\times n} vector (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs}). -%! -%! @strong{Remark 2.} If @code{dataset_.descriptive.cova} does not exist, the covariance matrix is computed prior to the -%! computation of the correlation matrix. -%! -%! @end deftypefn -%@eod: - -% Copyright (C) 2011-2012 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 . - -% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr - -if ~isfield(dataset_.descriptive,'cova') - dataset_ = compute_cova(dataset_); -end -normalization_matrix = diag(1./sqrt(diag(dataset_.descriptive.cova))); -dataset_.descriptive.corr = normalization_matrix*dataset_.descriptive.cova*normalization_matrix; \ No newline at end of file diff --git a/matlab/utilities/dataset/compute_cova.m b/matlab/utilities/dataset/compute_cova.m deleted file mode 100644 index 51f68981e..000000000 --- a/matlab/utilities/dataset/compute_cova.m +++ /dev/null @@ -1,67 +0,0 @@ -function dataset_ = compute_cova(dataset_) -% Computes the covariance matrix of the sample (possibly with missing observations). - -%@info: -%! @deftypefn {Function File} {@var{dataset_} =} compute_corr(@var{dataset_}) -%! @anchor{compute_corr} -%! This function computes covariance matrix of the sample (possibly with missing observations). -%! -%! @strong{Inputs} -%! @table @var -%! @item dataset_ -%! Dynare structure describing the dataset, built by @ref{initialize_dataset} -%! @end table -%! -%! @strong{Outputs} -%! @table @var -%! @item dataset_ -%! Dynare structure describing the dataset, built by @ref{initialize_dataset} -%! @end table -%! -%! @strong{This function is called by:} -%! @ref{descriptive_statistics}. -%! -%! @strong{This function calls:} -%! @ref{ndim}, @ref{demean}, @ref{nandemean}. -%! -%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.cova} is a -%! @tex{n\times n} vector (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs}). -%! -%! @end deftypefn -%@eod: - -% Copyright (C) 2011-2012 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 . - -% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr - -dataset_.descriptive.cova = zeros(dataset_.nvobs); - -data = transpose(dataset_.data); - -for i=1:dataset_.info.nvobs - for j=i:dataset_.info.nvobs - if dataset_.missing.state - dataset_.descriptive.cova(i,j) = nanmean(nandemean(data(:,i)).*nandemean(data(:,j))); - else - dataset_.descriptive.cova(i,j) = mean(demean(data(:,i)).*demean(data(:,j))); - end - if j>i - dataset_.descriptive.cova(j,i) = dataset_.descriptive.cova(i,j); - end - end -end \ No newline at end of file diff --git a/matlab/utilities/dataset/describe_missing_data.m b/matlab/utilities/dataset/describe_missing_data.m index 5506da2be..69589933e 100644 --- a/matlab/utilities/dataset/describe_missing_data.m +++ b/matlab/utilities/dataset/describe_missing_data.m @@ -26,9 +26,9 @@ function [i,n,s,j] = describe_missing_data(data) %! %! @end deftypefn %@eod: - -% Copyright (C) 2008-2012 Dynare Team -% + +% Copyright (C) 2008-2014 Dynare Team +% % This file is part of Dynare. % % Dynare is free software: you can redistribute it and/or modify @@ -85,7 +85,7 @@ end %@test:1 -%$ % Define a data set. +%$ % Define a data set. %$ A = [ 1 1 ; ... %$ 1 NaN ; ... %$ NaN 1 ; ... @@ -99,7 +99,7 @@ end %$ 1 1 ]; %$ %$ % Define expected results. -%$ eB = cell(1,11); +%$ eB = cell(1,11); %$ eB(1) = { transpose(1:2) }; %$ eB(2) = { 1 }; %$ eB(3) = { 2 }; @@ -107,7 +107,7 @@ end %$ eB(5) = { [] }; %$ eB(6) = { 1 }; %$ eB(7) = { 1 }; -%$ eB(8) = { transpose(1:2) }; +%$ eB(8) = { transpose(1:2) }; %$ eB(9) = { transpose(1:2) }; %$ eB(10) = { transpose(1:2) }; %$ eB(11) = { transpose(1:2) }; @@ -117,7 +117,7 @@ end %$ eE(1) = { [1; 2; 4; transpose(6:11)] }; %$ eE(2) = { [1; 3; 4; transpose(8:11)] }; %$ -%$ % Call the tested routine. +%$ % Call the tested routine. %$ [B,C,D,E] = describe_missing_data(transpose(A)); %$ %$ % Check the results. diff --git a/matlab/utilities/dataset/initialize_dataset.m b/matlab/utilities/dataset/initialize_dataset.m index 0b65cb8f3..5d992a288 100644 --- a/matlab/utilities/dataset/initialize_dataset.m +++ b/matlab/utilities/dataset/initialize_dataset.m @@ -1,4 +1,4 @@ -function dataset_ = initialize_dataset(datafile,varobs,first,nobs,transformation,prefilter,xls) +function dataset_ = initialize_dataset(datafile,varobs,first,nobs,logged_data_flag,prefilter,xls) % Initializes a structure describing the dataset. % Copyright (C) 2011-2013 Dynare Team @@ -18,8 +18,6 @@ function dataset_ = initialize_dataset(datafile,varobs,first,nobs,transformation % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr - if isempty(datafile) error('Estimation::initialize_dataset: You have to declare a dataset file!') end @@ -43,7 +41,7 @@ end % Fill the dataset structure dataset_.info.ntobs = nobs; -dataset_.info.nvobs = rows(varobs); +dataset_.info.nvobs = length(varobs); dataset_.info.varobs = varobs; % Test the number of variables in the database. @@ -61,9 +59,9 @@ if size(rawdata,1)~=dataset_.info.ntobs end rawdata = rawdata(first:(first+dataset_.info.ntobs-1),:); -% Take the log (or anything else) of the variables if needed -if isempty(transformation) - dataset_.rawdata = rawdata; +% Take the log of the variables if needed +if logged_data_flag + dataset_.rawdata = log(rawdata); else if isequal(transformation,@log) if ~isreal(rawdata) @@ -101,4 +99,4 @@ if prefilter == 1 dataset_.data = transpose(bsxfun(@minus,dataset_.rawdata,dataset_.descriptive.mean)); else dataset_.data = transpose(dataset_.rawdata); -end \ No newline at end of file +end diff --git a/matlab/utilities/dataset/makedataset.m b/matlab/utilities/dataset/makedataset.m new file mode 100644 index 000000000..9eb4cb396 --- /dev/null +++ b/matlab/utilities/dataset/makedataset.m @@ -0,0 +1,241 @@ +function [DynareDataset, DatasetInfo] = makedataset(DynareOptions, initialconditions, gsa_flag) + +% Initialize a dataset as a dseries object. +% +% +% INPUTS +% ====== +% +% DynareOptions [struct] Structure of options built by Dynare's preprocessor. +% +% +% OUTPUTS +% ======= +% +% DynareDataset [dseries] The dataset. +% DatasetInfo [struct] Various informations about the dataset (descriptive statistics and missing observations). +% +% EXAMPLE +% ======= +% +% [dataset_, dataset_info] = makedataset(options_) ; +% +% +% See also dynare_estimation_init + +if nargin<3 + gsa_flag = 0; +end + +if nargin<2 || isempty(initialconditions) + % If a the sample is to be used for the estimation of a VAR or DSGE-VAR model + % the second argument must be a strictly positive integer (the number of lags). + initialconditions = 0; +end + +if isempty(DynareOptions.datafile) && isempty(DynareOptions.dataset.file) && isempty(DynareOptions.dataset.series) + if gsa_flag + DynareDataset = dseries(); + DatasetInfo = []; + return + else + error('makedataset: datafile option is missing!') + end +end + +if isempty(DynareOptions.datafile) && ~isempty(DynareOptions.dataset.file) + datafile = DynareOptions.dataset.file; + newdatainterface = 1; +elseif isempty(DynareOptions.datafile) && ~isempty(DynareOptions.dataset.series) + try + dseriesobjectforuserdataset = evalin('base', DynareOptions.dataset.series); + catch + error(sprintf('makedataset: %s is unknown!', DynareOptions.dataset.series)) + end + if ~isdseries(dseriesobjectforuserdataset) + error(sprintf('makedataset: %s has to be a dseries object!', DynareOptions.dataset.series)) + end + datafile = []; + newdatainterface = 1; +elseif ~isempty(DynareOptions.datafile) && isempty(DynareOptions.dataset.file) + datafile = DynareOptions.datafile; + newdatainterface = 0; +elseif isempty(DynareOptions.datafile) && ~isempty(DynareOptions.dataset.file) + error('You cannot use simultaneously the data command and the datafile option (in the estimation command)!') +else + error('You have to specify the datafile!') +end + +% Check extension. +if ~isempty(datafile) + allowed_extensions = {'m','mat','csv','xls','xlsx'}; + datafile_extension = get_file_extension(datafile); + if isempty(datafile_extension) + available_extensions = {}; j = 1; + for i=1:length(allowed_extensions) + if exist([datafile '.' allowed_extensions{i}]) + available_extensions(j) = {allowed_extensions{i}}; + j = j+1; + end + end + if isempty(available_extensions) + error(['I can''t find a datafile (with allowed extension)!']) + end + if length(available_extensions)>1 + error(sprintf(['You did not specify an extension for the datafile, but more than one candidate ' ... + 'are available in the designed folder!\nPlease, add an extension to the datafile ' ... + '(m, mat, csv, xls or xlsx are legal extensions).'])); + end + datafile = [datafile '.' available_extensions{1}]; + end +end + +% Load the data in a dseries object. +if ~isempty(datafile) + DynareDataset = dseries(datafile); +else + DynareDataset = dseriesobjectforuserdataset; + clear('dseriesobjectforuserdataset'); +end + +% Select a subset of the variables. +DynareDataset = DynareDataset{DynareOptions.varobs{:}}; + +% Apply log function if needed. +if DynareOptions.loglinear && ~DynareOptions.logdata + DynareDataset = DynareDataset.log(); +end + +% Test if an initial period (different from its default value) is explicitely defined in the datafile. +if isequal(DynareDataset.init, dates(1,1)) + dataset_default_initial_period = 1; +else + dataset_default_initial_period = 0; +end + +% Test if an initial period (different from its default value) is explicitely defined in the mod file with the set_time command. +if isequal(DynareOptions.initial_period, dates(1,1)) + set_time_default_initial_period = 1; +else + set_time_default_initial_period = 0; +end + +if ~set_time_default_initial_period && dataset_default_initial_period + % Overwrite the initial period in dataset (it was set to default). + % Note that the updates of freq and time members are auto-magically + % done by dseries::subsasgn overloaded method. + DynareDataset.init = DynareOptions.initial_period; +end + +if set_time_default_initial_period && ~dataset_default_initial_period + % Overwrite the global initial period defined by set_time (it was set to default). + DynareOptions.initial_period = DynareDataset.init; +end + +if ~set_time_default_initial_period && ~dataset_default_initial_period + % Check if dataset.init and options_.initial_period are identical. + if DynareOptions.initial_period>DynareDataset.init + %if ~isequal(DynareOptions.initial_period, DynareDataset.init) + error('The date as defined by the set_time command is not consistent with the initial period in the database!') + end +end + +% Set firstobs, lastobs and nobs +if newdatainterface + if isempty(DynareOptions.dataset.firstobs) + % first_obs option was not used in the data command. + firstobs = DynareDataset.init; + else + firstobs = DynareOptions.dataset.firstobs; + end + if isnan(DynareOptions.dataset.nobs) + % nobs option was not used in the data command. + if isempty(DynareOptions.dataset.lastobs) + % last_obs option was not used in the data command. + nobs = DynareDataset.nobs; + lastobs = DynareDataset.dates(end); + else + lastobs = DynareOptions.dataset.lastobs; + nobs = lastobs-firstobs+1; + end + else + nobs = DynareOptions.dataset.nobs; + if isempty(DynareOptions.dataset.lastobs) + % last_obs option was not used in the data command. + lastobs = firstobs+(nobs-1); + else + % last_obs and nobs were used in the data command. Check that they are consistent (with firstobs). + if ~isequal(lastobs,firstobs+(nobs-1)) + error(sprintf('Options last_obs (%s), first_obs (%s) and nobs (%s) are not consistent!',char(lastobs),char(firstobs),num2str(nobs))); + end + end + end +else + if isnan(DynareOptions.first_obs) + firstobs = DynareDataset.init; + else + firstobs = DynareDataset.dates(DynareOptions.first_obs); + end + if isnan(DynareOptions.nobs) + lastobs = DynareDataset.dates(end); + nobs = lastobs-firstobs+1; + else + nobs = DynareOptions.nobs; + lastobs = firstobs+(nobs-1); + end +end + +% Add initial conditions if needed +FIRSTOBS = firstobs-initialconditions; + +% Check that firstobs belongs to DynareDataset.dates +if firstobsDynareDataset.dates(end) + error(sprintf('last_obs (%s) cannot be greater than the last date in the dataset (%s)!',char(lastobs),char(DynareDataset.dates(end)))) + end +else + % ... or check that nobs is smaller than the number of observations in DynareDataset. + if nobs>DynareDataset.nobs + error(sprintf('nobs (%s) cannot be greater than the last date in the dataset (%s)!', num2str(nobs), num2str(nobs))) + end +end + +% Select a subsample. +DynareDataset = DynareDataset(FIRSTOBS:lastobs); + +% Initialize DatasetInfo structure. +DatasetInfo = struct('missing', struct('state', NaN, 'aindex', [], 'vindex', [], 'number_of_observations', NaN, 'no_more_missing_observations', NaN), ... + 'descriptive', struct('mean', [], 'covariance', [], 'correlation', [], 'autocovariance', [])); + +% Fill DatasetInfo.missing if some observations are missing +DatasetInfo.missing.state = isanynan(DynareDataset.data); +if DatasetInfo.missing.state + [DatasetInfo.missing.aindex, DatasetInfo.missing.number_of_observations, DatasetInfo.missing.no_more_missing_observations, DatasetInfo.missing.vindex] = ... + describe_missing_data(DynareDataset.data); +else + DatasetInfo.missing.aindex = num2cell(transpose(repmat(1:DynareDataset.vobs,DynareDataset.nobs,1)),1); + DatasetInfo.missing.no_more_missing_observations = 1; +end + +% Compute the empirical mean of the observed variables. +DatasetInfo.descriptive.mean = nanmean(DynareDataset.data); + +% Compute the empirical covariance matrix of the observed variables. +DatasetInfo.descriptive.covariance = nancovariance(DynareDataset.data); + +% Compute the empirical correlation matrix of the observed variables. +normalization_matrix = diag(1./sqrt(diag(DatasetInfo.descriptive.covariance))); +DatasetInfo.descriptive.correlation = normalization_matrix*DatasetInfo.descriptive.covariance*normalization_matrix; + +% Compute autocorrelation function. +DatasetInfo.descriptive.autocovariance = nanautocovariance(DynareDataset.data, DynareOptions.ar); \ No newline at end of file diff --git a/matlab/utilities/dataset/compute_acov.m b/matlab/utilities/dataset/nanautocovariance.m similarity index 57% rename from matlab/utilities/dataset/compute_acov.m rename to matlab/utilities/dataset/nanautocovariance.m index b1834e854..715ae957f 100644 --- a/matlab/utilities/dataset/compute_acov.m +++ b/matlab/utilities/dataset/nanautocovariance.m @@ -1,4 +1,4 @@ -function dataset_ = compute_acov(dataset_) +function autocov = nanautocovariance(data,order) % Computes the (multivariate) auto-covariance function of the sample (possibly with missing observations). %@info: @@ -8,36 +8,36 @@ function dataset_ = compute_acov(dataset_) %! %! @strong{Inputs} %! @table @var -%! @item dataset_ -%! Dynare structure describing the dataset, built by @ref{initialize_dataset} -%! @item nlag -%! Integer scalar. The maximum number of lags of the autocovariance function. +%! @item data +%! T*N array of real numbers. +%! @item order +%! Integer scalar. The maximum number of lags of the autocovariance function. %! @end table %! %! @strong{Outputs} %! @table @var -%! @item dataset_ -%! Dynare structure describing the dataset, built by @ref{initialize_dataset} +%! @item autocov +%! A N*N*order array of real numbers. %! @end table -%! -%! @strong{This function is called by:} +%! +%! @strong{This function is called by:} %! @ref{descriptive_statistics}. -%! +%! %! @strong{This function calls:} -%! @ref{ndim}, @ref{compute_cova}, @ref{demean}, @ref{nandemean}. -%! -%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.acov} is a -%! @tex{n\times n\times p} array (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs}, +%! @ref{ndim}, @ref{nancovariance}, @ref{demean}, @ref{nandemean}. +%! +%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.acov} is a +%! @tex{n\times n\times p} array (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs}, %! and @tex{n} is the maximum number of lags given by the second input @code{nlag}). -%! -%! @strong{Remark 2.} If @code{dataset_.descriptive.cova} does not exist, the covariance matrix is computed prior to the +%! +%! @strong{Remark 2.} If @code{dataset_.descriptive.cova} does not exist, the covariance matrix is computed prior to the %! computation of the auto-covariance function. %! %! @end deftypefn %@eod: -% Copyright (C) 2011-2012 Dynare Team -% +% Copyright (C) 2011-2014 Dynare Team +% % This file is part of Dynare. % % Dynare is free software: you can redistribute it and/or modify @@ -53,22 +53,23 @@ function dataset_ = compute_acov(dataset_) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr +n = size(data,2); +missing = isanynan(data); -if ~isfield(dataset_.descriptive,'cova') - dataset_ = compute_cova(dataset_); -end -dataset_.descriptive.acov = zeros(dataset_.nvobs,dataset_.nvobs,nlag); +autocov = zeros(n, n, order); -data = transpose(dataset_.data); - -for lag=1:nlag - for i=1:dataset_.info.nvobs - for j=1:dataset_.info.nvobs - if dataset_.missing.state - dataset_.descriptive.acov(i,j,lag) = nanmean(nandemean(data(lag+1:end,i)).*nandemean(data(1:end-lag,j))); +for lag=1:order + if missing + data = nandemean(data); + else + data = demean(data); + end + for i=1:n + for j=1:n + if missing + autocov(i,j,lag) = nanmean(data((lag+1):end,i).*data(1:end-lag,j)); else - dataset_.descriptive.acov(i,j,lag) = mean(demean(data(lag+1:end,i)).*demean(data(1:end-lag,j))); + autocov(i,j,lag) = mean(data((lag+1):end,i).*data(1:end-lag,j)); end end end diff --git a/matlab/utilities/dataset/nancovariance.m b/matlab/utilities/dataset/nancovariance.m new file mode 100644 index 000000000..8e2e2b227 --- /dev/null +++ b/matlab/utilities/dataset/nancovariance.m @@ -0,0 +1,101 @@ +function CovarianceMatrix = nancovariance(data) +% Computes the covariance matrix of a sample (possibly with missing observations). + +%@info: +%! @deftypefn {Function File} {@var{CovarianceMatrix} =} compute_corr(@var{data}) +%! @anchor{compute_cova} +%! This function computes covariance matrix of a sample defined by a dseries object (possibly with missing observations). +%! +%! @strong{Inputs} +%! @table @var +%! @item data +%! a T*N array of real numbers. +%! @end table +%! +%! @strong{Outputs} +%! @table @var +%! @item CovarianceMatrix +%! Array of real numbers. +%! @end table +%! +%! @strong{This function is called by:} +%! @ref{descriptive_statistics}. +%! +%! @strong{This function calls:} +%! @ref{ndim}, @ref{demean}, @ref{nandemean}. +%! +%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.cova} is a +%! @tex{n\times n} vector (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs}). +%! +%! @end deftypefn +%@eod: + +% Copyright (C) 2011-2014 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 . + +% Initialize the output. +CovarianceMatrix = zeros(size(data,2)); + +if isanynan(data) + data = bsxfun(@minus,data,nanmean(data)); + for i=1:size(data,2) + for j=i:size(data,2) + CovarianceMatrix(i,j) = nanmean(data(:,i).*data(:,j)); + if j>i + CovarianceMatrix(j,i) = CovarianceMatrix(i,j); + end + end + end +else + data = bsxfun(@minus,data,mean(data)); + CovarianceMatrix = (transpose(data)*data)/size(data,1); +end + +%@test:1 +%$ +%$ % Define a dataset. +%$ data1 = randn(10000000,2); +%$ +%$ % Same dataset with missing observations. +%$ data2 = data1; +%$ data2(45,1) = NaN; +%$ data2(57,2) = NaN; +%$ data2(367,:) = NaN(1,2); +%$ +%$ t = zeros(2,1); +%$ +%$ % Call the tested routine. +%$ try +%$ c1 = nancovariance(data1); +%$ t(1) = 1; +%$ catch +%$ t(1) = 0; +%$ end +%$ try +%$ c2 = nancovariance(data2); +%$ t(2) = 1; +%$ catch +%$ t(2) = 0; +%$ end +%$ +%$ if t(1) && t(2) +%$ t(3) = max(max(abs(c1-c2)))<1e-4; +%$ end +%$ +%$ % Check the results. +%$ T = all(t); +%@eof:1 \ No newline at end of file diff --git a/matlab/utilities/dataset/nanmoments.m b/matlab/utilities/dataset/nanmoments.m new file mode 100644 index 000000000..78e95e9b0 --- /dev/null +++ b/matlab/utilities/dataset/nanmoments.m @@ -0,0 +1,25 @@ +function m = nanmoments(data, n) +% Compute centered marginal moments of order n (possibly with missing observations). + +% Copyright (C) 2014 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 . + +if isanynan(data) + m = transpose(nanmean(bsxfun(@power,nandemean(data),n))); +else + m = transpose(mean(bsxfun(@power,demean(data),n))); +end \ No newline at end of file diff --git a/matlab/utilities/dataset/compute_stdv.m b/matlab/utilities/dataset/nanvariance.m similarity index 62% rename from matlab/utilities/dataset/compute_stdv.m rename to matlab/utilities/dataset/nanvariance.m index 4384f83dd..d4894e8e9 100644 --- a/matlab/utilities/dataset/compute_stdv.m +++ b/matlab/utilities/dataset/nanvariance.m @@ -1,21 +1,21 @@ -function dataset_ = compute_stdv(dataset_) +function variances = nanvariance(data) % Compute the standard deviation for each observed variable (possibly with missing observations). %@info: -%! @deftypefn {Function File} {@var{dataset_} =} compute_stdv(@var{dataset_}) -%! @anchor{compute_stdv} -%! This function computes the standard deviation of the observed variables (possibly with missing observations). +%! @deftypefn {Function File} {@var{variances} =} nanvariance(@var{data}) +%! @anchor{nanvariance} +%! This function computes the variances of the observed variables (possibly with missing observations). %! %! @strong{Inputs} %! @table @var -%! @item dataset_ -%! Dynare structure describing the dataset, built by @ref{initialize_dataset} +%! @item datas +%! A T*N array of real numbers. %! @end table %! %! @strong{Outputs} %! @table @var -%! @item dataset_ -%! Dynare structure describing the dataset, built by @ref{initialize_dataset} +%! @item variances +%! A N*1 vector of real numbers %! @end table %! %! @strong{This function is called by:} @@ -30,7 +30,7 @@ function dataset_ = compute_stdv(dataset_) %! @end deftypefn %@eod: -% Copyright (C) 2011-2012 Dynare Team +% Copyright (C) 2011-2014 Dynare Team % % This file is part of Dynare. % @@ -47,10 +47,8 @@ function dataset_ = compute_stdv(dataset_) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr - -if dataset_.missing.state - dataset_.descriptive.stdv = sqrt(nanmean(bsxfun(@power,nandemean(transpose(dataset_.data)),2))); +if isanynan(data) + variances = transpose(nanmean(bsxfun(@power,nandemean(data),2))); else - dataset_.descriptive.stdv = sqrt(mean(bsxfun(@power,demean(transpose(dataset_.data)),2))); + variances = transpose(mean(bsxfun(@power,demean(data),2))); end \ No newline at end of file diff --git a/matlab/utilities/dates/isdate.m b/matlab/utilities/dates/isdate.m index 2fb603803..726e5e479 100644 --- a/matlab/utilities/dates/isdate.m +++ b/matlab/utilities/dates/isdate.m @@ -25,15 +25,16 @@ function b = isdate(str) % --*-- Unitary tests --*-- % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -%if length(str)>1 - b = isquaterly(str) || isyearly(str) || ismonthly(str) || isweekly(str); - %else - %b = 0; - %end +if isnumeric(str) && isscalar(str) + b = 1; + return +end + +b = isstringdate(str); %@test:1 %$ -%$ date_1 = '1950M2'; +%$ date_1 = 1950; %$ date_2 = '1950m2'; %$ date_3 = '-1950m2'; %$ date_4 = '1950m52'; diff --git a/matlab/utilities/dates/isfreq.m b/matlab/utilities/dates/isfreq.m index 38a779f20..79ba42c53 100644 --- a/matlab/utilities/dates/isfreq.m +++ b/matlab/utilities/dates/isfreq.m @@ -25,4 +25,15 @@ function B = isfreq(A) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -B = (ischar(A) && ismember(upper(A),{'Y','A','Q','M','W'})) || ismember(A,[1 4 12 52]); \ No newline at end of file +B = 0; + +if ischar(A) + if isequal(length(A),1) && ismember(upper(A),{'Y','A','Q','M','W'}) + B = 1; + return + end +end + +if isnumeric(A) && isequal(length(A),1) && ismember(A,[1 4 12 52]) + B = 1; +end \ No newline at end of file diff --git a/matlab/utilities/dates/isstringdate.m b/matlab/utilities/dates/isstringdate.m new file mode 100644 index 000000000..fbac8b007 --- /dev/null +++ b/matlab/utilities/dates/isstringdate.m @@ -0,0 +1,54 @@ +function b = isstringdate(str) % --*-- Unitary tests --*-- + +% Tests if the input string can be interpreted as a date. +% +% INPUTS +% o str string. +% +% OUTPUTS +% o b integer scalar, equal to 1 if str can be interpreted as a date or 0 otherwise. + +% Copyright (C) 2013 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 . + +if ischar(str) + b = isquaterly(str) || isyearly(str) || ismonthly(str) || isweekly(str); +else + b = 0; +end + +%@test:1 +%$ +%$ date_1 = '1950M2'; +%$ date_2 = '1950m2'; +%$ date_3 = '-1950m2'; +%$ date_4 = '1950m52'; +%$ date_5 = ' 1950'; +%$ date_6 = '1950Y'; +%$ date_7 = '-1950a'; +%$ date_8 = '1950m '; +%$ +%$ t(1) = dyn_assert(isstringdate(date_1),1); +%$ t(2) = dyn_assert(isstringdate(date_2),1); +%$ t(3) = dyn_assert(isstringdate(date_3),1); +%$ t(4) = dyn_assert(isstringdate(date_4),0); +%$ t(5) = dyn_assert(isstringdate(date_5),0); +%$ t(6) = dyn_assert(isstringdate(date_6),1); +%$ t(7) = dyn_assert(isstringdate(date_7),1); +%$ t(8) = dyn_assert(isstringdate(date_8),0); +%$ T = all(t); +%@eof:1 \ No newline at end of file diff --git a/matlab/utilities/dseries/from.m b/matlab/utilities/dseries/from.m index a8ae4bef3..02be01367 100644 --- a/matlab/utilities/dseries/from.m +++ b/matlab/utilities/dseries/from.m @@ -1,4 +1,4 @@ -function from(varargin) +function from(varargin) % --*-- Unitary tests --*-- % Copyright (C) 2014 Dynare Team % @@ -219,10 +219,10 @@ for i=1:number_of_variables end if i>1 if ismember(var.name,variable_names) - error('dseries::from: All the dseries objects should contain variables with different names!') - else - variable_names(i) = {var.name{1}}; + % Locally change variable name. + var = var.rename(var.name{1},get_random_string(20)); end + variable_names(i) = {var.name{1}}; else variable_names(i) = {var.name{1}}; end @@ -421,4 +421,18 @@ function i = isassignedvariable(var,expr) return end end - i = 0; \ No newline at end of file + i = 0; + +%@test:1 +%$ try +%$ y = dseries(zeros(400,1),dates('1950Q1')) ; +%$ v = dseries(randn(400,1),dates('1950Q1')) ; +%$ u = dseries(randn(400,1),dates('1950Q1')) ; +%$ from 1950Q2 to 2049Q4 do y(t) = (1+.01*u(t))*y(t-1) + v(t) +%$ t(1) = 1; +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ T = all(t); +%@eof:1 diff --git a/matlab/utilities/general/isanynan.m b/matlab/utilities/general/isanynan.m new file mode 100644 index 000000000..7a16419d0 --- /dev/null +++ b/matlab/utilities/general/isanynan.m @@ -0,0 +1,21 @@ +function yes = isanynan(array) +% Return one if the array contains at least one NaN, 0 otherwise. + +% Copyright (C) 2011-2014 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 . + +yes = any(isnan(array(:))); \ No newline at end of file diff --git a/matlab/var_sample_moments.m b/matlab/var_sample_moments.m index ab7b4b74b..f9d4f8b02 100644 --- a/matlab/var_sample_moments.m +++ b/matlab/var_sample_moments.m @@ -1,29 +1,28 @@ -function [YtY,XtY,YtX,XtX,Y,X] = ... - var_sample_moments(FirstObservation,LastObservation,qlag,var_trend_order,datafile,varobs,xls_sheet,xls_range) +function var_sample_moments(nlag, var_trend_order, dataset_)%datafile,varobs,xls_sheet,xls_range) % Computes the sample moments of a VAR model. % % The VAR(p) model is defined by: % -% y_t = \sum_{k=1}^p y_{t-k} A_k + z_t C + e_t for t = 1,...,T +% y_t = \sum_{k=1}^p y_{t-k} A_k + z_t C + e_t for t = 1,...,T % % where y_t is a 1*m vector of observed endogenous variables, p is the % number of lags, A_k is an m*m real matrix, z_t is a 1*q vector of % exogenous (deterministic) variables, C is a q*m real matrix and % e_t is a vector of exogenous stochastic shocks. T is the number -% of observations. The deterministic exogenous variables are assumed to -% be a polynomial trend of order q = "var_trend_order". +% of observations. The deterministic exogenous variables are assumed to +% be a polynomial trend of order q = "var_trend_order". % -% We define: +% We define: % % <> Y = (y_1',y_2',...,y_T')' a T*m matrix, % -% <> x_t = (y_{t-1},y_{t-2},...,y_{t-p},z_t) a 1*(mp+q) row vector, +% <> x_t = (y_{t-1},y_{t-2},...,y_{t-p},z_t) a 1*(mp+q) row vector, % -% <> X = (x_1',x_2',...,x_T')' a T*(mp+q) matrix, +% <> X = (x_1',x_2',...,x_T')' a T*(mp+q) matrix, % % <> E = (e_1',e_2',...,e_T')' a T*m matrix and % -% <> A = (A_1',A_2',...,A_p',C')' an (mp+q)*m matrix of coefficients. +% <> A = (A_1',A_2',...,A_p',C')' an (mp+q)*m matrix of coefficients. % % So that we can equivalently write the VAR(p) model using the following % matrix representation: @@ -31,18 +30,17 @@ function [YtY,XtY,YtX,XtX,Y,X] = ... % Y = X * A +E % % -% INPUTS -% o FirstObservation [integer] First observation. -% o LastObservation [integer] Last observation. -% o qlag [integer] Number of lags in the VAR model. -% o var_trend_order [integer] Order of the polynomial exogenous trend: +% INPUTS +% o nlag [integer] Number of lags in the VAR model. +% o var_trend_order [integer] Order of the polynomial exogenous trend: % = -1 no constant and no linear trend, % = 0 constant and no linear trend, % = 1 constant and linear trend. +% o dataset_ [dseries] The sample. % -% OUTPUTS +% OUTPUTS % o YtY [double] Y'*Y an m*m matrix. -% o XtY [double] X'*Y an (mp+q)*m matrix. +% o XtY [double] X'*Y an (mp+q)*m matrix. % o YtX [double] Y'*X an m*(mp+q) matrix. % o XtX [double] X'*X an (mp+q)*(mp+q) matrix. % o Y [double] Y a T*m matrix. @@ -50,8 +48,11 @@ function [YtY,XtY,YtX,XtX,Y,X] = ... % % SPECIAL REQUIREMENTS % None. +% +% REMARKS +% Outputs are saved in the base workspace (not returned by the function). -% Copyright (C) 2007-2009 Dynare Team +% Copyright (C) 2007-2014 Dynare Team % % This file is part of Dynare. % @@ -68,42 +69,35 @@ function [YtY,XtY,YtX,XtX,Y,X] = ... % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -X = []; -Y = []; -YtY = []; -YtX = []; -XtY = []; -XtX = []; +LastObservation = dataset_.dates(end); +FirstObservation = dataset_.dates(1)+nlag; -data = read_variables(datafile,varobs,[],xls_sheet,xls_range); +NumberOfObservations = LastObservation-FirstObservation+1; +NumberOfVariables = dataset_.vobs; -if qlag > FirstObservation - error('VarSampleMoments :: not enough data to initialize! Try to increase FirstObservation.') - return -end - -NumberOfObservations = LastObservation-FirstObservation+1;% This is T. -NumberOfVariables = size(varobs,1);% This is m. -if var_trend_order == -1% No constant no linear trend case. - X = zeros(NumberOfObservations,NumberOfVariables*qlag); -elseif var_trend_order == 0% Constant and no linear trend case. -X = ones(NumberOfObservations,NumberOfVariables*qlag+1); -indx = NumberOfVariables*qlag+1; -elseif var_trend_order == 1;% Constant and linear trend case. -X = ones(NumberOfObservations,NumberOfVariables*qlag+2); -indx = NumberOfVariables*qlag+1:NumberOfVariables*qlag+2; +if isequal(var_trend_order,-1) + % No constant no linear trend case. + X = zeros(NumberOfObservations,NumberOfVariables*nlag); +elseif isequal(var_trend_order,0) + % Constant and no linear trend case. + X = ones(NumberOfObservations,NumberOfVariables*nlag+1); + indx = NumberOfVariables*nlag+1; +elseif isequal(var_trend_order,1) + % Constant and linear trend case. + X = ones(NumberOfObservations,NumberOfVariables*nlag+2); + indx = NumberOfVariables*nlag+1:NumberOfVariables*nlag+2; else - error('var_sample_moments :: trend must be equal to -1,0 or 1!') + error('Estimation::var_sample_moments: trend must be equal to -1,0 or 1!') return end -% I build matrices Y and X -Y = data(FirstObservation:LastObservation,:); +% I build matrices Y and X +Y = dataset_(FirstObservation:LastObservation).data; for t=1:NumberOfObservations - line = t + FirstObservation-1; - for lag = 1:qlag - X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = data(line-lag,:); + currentdate = FirstObservation+(t-1); + for lag = 1:nlag + X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = dataset_(currentdate-lag).data; end end @@ -111,7 +105,9 @@ if (var_trend_order == 1) X(:,end) = transpose(1:NumberOfObservations) end -YtY = Y'*Y; -YtX = Y'*X; -XtY = X'*Y; -XtX = X'*X; \ No newline at end of file +assignin('base', 'mYY', Y'*Y); +assignin('base', 'mYX', Y'*X); +assignin('base', 'mXY', X'*Y); +assignin('base', 'mXX', X'*X); +assignin('base', 'Ydata', Y); +assignin('base', 'Xdata', X); \ No newline at end of file diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index e9292f901..c23908a0d 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -1876,9 +1876,17 @@ EstimationDataStatement::checkPass(ModFileStructure &mod_file_struct, WarningCon exit(EXIT_FAILURE); } - if (options_list.string_options.find("file") == options_list.string_options.end()) + if ((options_list.string_options.find("file") == options_list.string_options.end()) && + (options_list.string_options.find("series") == options_list.string_options.end())) { - cerr << "ERROR: The file option must be passed to the data statement." << endl; + cerr << "ERROR: The file or series option must be passed to the data statement." << endl; + exit(EXIT_FAILURE); + } + + if ((options_list.string_options.find("file") != options_list.string_options.end()) && + (options_list.string_options.find("series") != options_list.string_options.end())) + { + cerr << "ERROR: The file and series options cannot be used simultaneously in the data statement." << endl; exit(EXIT_FAILURE); } } @@ -1887,8 +1895,8 @@ void EstimationDataStatement::writeOutput(ostream &output, const string &basename) const { options_list.writeOutput(output, "options_.dataset"); - if (options_list.date_options.find("first_obs") == options_list.date_options.end()) - output << "options_.dataset.firstobs = options_.initial_period;" << endl; + //if (options_list.date_options.find("first_obs") == options_list.date_options.end()) + // output << "options_.dataset.first_obs = options_.initial_period;" << endl; } SubsamplesStatement::SubsamplesStatement(const string &name1_arg, diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index c2cb723e7..0c754a4bc 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -90,7 +90,7 @@ class ParsingDriver; %token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED %token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION %token CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED -%token DATAFILE FILE DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS +%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR %token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME %token FLOAT_NUMBER DATES @@ -833,7 +833,7 @@ restriction_expression : expression {driver.check_restriction_expression_constan restriction_expression_1 : restriction_elem_expression | restriction_expression_1 restriction_elem_expression - ; + ; restriction_elem_expression : COEFF '(' symbol COMMA INT_NUMBER ')' { driver.add_positive_restriction_element($3,$5);} @@ -1366,9 +1366,10 @@ data_options_list : data_options_list COMMA data_options ; data_options : o_file - | o_new_estimation_data_first_obs - | o_last_obs - | o_new_estimation_data_nobs + | o_series + | o_data_first_obs + | o_data_last_obs + | o_data_nobs | o_xls_sheet | o_xls_range ; @@ -2493,6 +2494,7 @@ o_simul_seed : SIMUL_SEED EQUAL INT_NUMBER { driver.error("'simul_seed' option i o_qz_criterium : QZ_CRITERIUM EQUAL non_negative_number { driver.option_num("qz_criterium", $3); }; o_qz_zero_threshold : QZ_ZERO_THRESHOLD EQUAL non_negative_number { driver.option_num("qz_zero_threshold", $3); }; o_file : FILE EQUAL filename { driver.option_str("file", $3); }; +o_series : SERIES EQUAL symbol { driver.option_str("series", $3); }; o_datafile : DATAFILE EQUAL filename { driver.option_str("datafile", $3); }; o_nobs : NOBS EQUAL vec_int { driver.option_vec_int("nobs", $3); } @@ -2505,12 +2507,9 @@ o_conditional_variance_decomposition : CONDITIONAL_VARIANCE_DECOMPOSITION EQUAL { driver.option_vec_int("conditional_variance_decomposition", $3); } ; o_first_obs : FIRST_OBS EQUAL INT_NUMBER { driver.option_num("first_obs", $3); }; -o_new_estimation_data_first_obs : FIRST_OBS EQUAL date_expr - { driver.option_date("first_obs", $3); } - ; -o_last_obs : LAST_OBS EQUAL date_expr - { driver.option_date("last_obs", $3); } - ; +o_data_first_obs : FIRST_OBS EQUAL date_expr { driver.option_date("firstobs", $3); } ; +o_data_last_obs : LAST_OBS EQUAL date_expr { driver.option_date("lastobs", $3); } ; +o_data_nobs : NOBS EQUAL INT_NUMBER { driver.option_num("nobs", $3); }; o_shift : SHIFT EQUAL signed_number { driver.option_num("shift", $3); }; o_shape : SHAPE EQUAL prior_distribution { driver.prior_shape = $3; }; o_mode : MODE EQUAL signed_number { driver.option_num("mode", $3); }; @@ -2523,7 +2522,6 @@ o_bounds : BOUNDS EQUAL vec_value_w_inf { driver.option_num("bounds", $3); }; o_domain : DOMAINN EQUAL vec_value { driver.option_num("domain", $3); }; o_interval : INTERVAL EQUAL vec_value { driver.option_num("interval", $3); }; o_variance : VARIANCE EQUAL expression { driver.set_prior_variance($3); } -o_new_estimation_data_nobs : NOBS EQUAL INT_NUMBER { driver.option_num("nobs", $3); }; o_prefilter : PREFILTER EQUAL INT_NUMBER { driver.option_num("prefilter", $3); }; o_presample : PRESAMPLE EQUAL INT_NUMBER { driver.option_num("presample", $3); }; o_lik_algo : LIK_ALGO EQUAL INT_NUMBER { driver.option_num("lik_algo", $3); }; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 6ab6faab0..a552f1997 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -604,6 +604,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2 simul_replic {return token::SIMUL_REPLIC;} xls_sheet {return token::XLS_SHEET;} xls_range {return token::XLS_RANGE;} +series {return token::SERIES;} mh_recover {return token::MH_RECOVER;} planner_discount {return token::PLANNER_DISCOUNT;} calibration {return token::CALIBRATION;} diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 70889c6b2..db02db939 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -135,7 +135,20 @@ OptionsList::writeOutput(ostream &output) const void OptionsList::writeOutput(ostream &output, const string &option_group) const { - output << option_group << " = struct();" << endl; + // Initialize option_group as an empty struct iff the field does not exist! + unsigned idx = option_group.find_last_of("."); + if (idx 0) { - output << "options_.varobs = [];" << endl; + int ic = 1; + output << "options_.varobs = cell(1);" << endl; for (vector::const_iterator it = varobs.begin(); it != varobs.end(); it++) - if (it == varobs.begin()) - output << "options_.varobs = '" << getName(*it) << "';" << endl; - else - output << "options_.varobs = char(options_.varobs, '" << getName(*it) << "');" << endl; - + { + output << "options_.varobs(" << ic << ") = {'" << getName(*it) << "'};" << endl; + ic++; + } output << "options_.varobs_id = [ "; for (vector::const_iterator it = varobs.begin(); it != varobs.end(); it++) diff --git a/tests/Makefile.am b/tests/Makefile.am index 2e8b02fa9..60ac51c55 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -80,6 +80,7 @@ MODFILES = \ fs2000/fs2000e.mod \ fs2000/fs2000_cmaes.mod \ fs2000/fs2000_calib.mod \ + fs2000/fs2000_calib_dseries.mod \ fs2000/fs2000_analytic_derivation.mod \ fs2000/fs2000_missing_data.mod \ fs2000/fs2000_sd.mod \ diff --git a/tests/fs2000/fs2000_calib_dseries.mod b/tests/fs2000/fs2000_calib_dseries.mod new file mode 100644 index 000000000..4354acf7e --- /dev/null +++ b/tests/fs2000/fs2000_calib_dseries.mod @@ -0,0 +1,73 @@ +// See fs2000.mod in the examples/ directory for details on the model + +set_time(1950Q1); + +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; + + +varobs gp_obs gy_obs; + +fsdata = dseries('fsdat_simul.m',1950Q1); +fsdata = fsdata{'g@p,y@_obs'}(1960Q1:1980Q4); + +fsdata + +data(series=fsdata);//, first_obs=1960Q1, last_obs=1980Q4); + + +calib_smoother(filtered_vars, filter_step_ahead = [3:4]) m P c e W R k d n l y dA; + diff --git a/tests/fs2000/fs2000_data.mod b/tests/fs2000/fs2000_data.mod new file mode 100644 index 000000000..23018eedd --- /dev/null +++ b/tests/fs2000/fs2000_data.mod @@ -0,0 +1,79 @@ +// See fs2000.mod in the examples/ directory for details on the model + +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; + +options_.solve_tolf = 1e-12; + +set_time(1970Q3); // Interpreted as the first date available in the sample loaded below. + +data(file='fsdat_simul.m',first_obs=1971Q1, nobs=40); + +estimation(order=1,nobs=192,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8); diff --git a/tests/fs2000/fs2000_dseries_a.mod b/tests/fs2000/fs2000_dseries_a.mod new file mode 100644 index 000000000..1210e7425 --- /dev/null +++ b/tests/fs2000/fs2000_dseries_a.mod @@ -0,0 +1,77 @@ +set_time(1950Q3); + +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; + +options_.solve_tolf = 1e-12; + +data(file=fsdat_simul_dseries,first_obs=1950Q3, nobs=192); + +estimation(order=1,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8); diff --git a/tests/fs2000/fs2000_dseries_b.mod b/tests/fs2000/fs2000_dseries_b.mod new file mode 100644 index 000000000..c2e31087c --- /dev/null +++ b/tests/fs2000/fs2000_dseries_b.mod @@ -0,0 +1,80 @@ +set_time(1950Q3); + +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; + +options_.solve_tolf = 1e-12; + +fsdataset = dseries('fsdat_simul_dseries.m'); +fsdataset = fsdataset(1950Q3:1950Q3+191); + +data(series=fsdataset); + +estimation(order=1,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8); diff --git a/tests/gsa/ls2003.mod b/tests/gsa/ls2003.mod index 883eab9ef..6cf6410b0 100644 --- a/tests/gsa/ls2003.mod +++ b/tests/gsa/ls2003.mod @@ -108,7 +108,7 @@ disp(' '); disp('MC FILTERING(rmse=1), TO MAP THE FIT FROM PRIORS'); disp('Press ENTER to continue'); pause(5); -dynare_sensitivity(nodisplay, datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, // also presample=2,loglinear, are admissible +dynare_sensitivity(nodisplay, datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1, // also presample=2,loglinear, are admissible load_stab=1, // load prior sample istart_rmse=2, //start computing rmse from second observation (i.e. rmse does not inlude initial big error) stab=0, // don't plot again stability analysis results @@ -121,7 +121,7 @@ disp('BY USING THE COMBINED CALL'); disp(' '); disp('dynare_sensitivity(redform=1,') disp('logtrans_redform=1, namendo=(pie,R), namexo=(e_R), namlagendo=(R),') -disp('datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,') +disp('datafile=data_ca1.m,first_obs=8,nobs=79,prefilter=1,') disp('istart_rmse=2, rmse=1);') disp(' '); disp('Press ENTER to continue'); pause(5); @@ -131,7 +131,7 @@ disp('Press ENTER to continue'); pause(5); //namendo=(pie,R), // evaluate relationships for pie and R (namendo=(:) for all variables) //namexo=(e_R), // evaluate relationships with exogenous e_R (use namexo=(:) for all shocks) //namlagendo=(R), // evaluate relationships with lagged R (use namlagendo=(:) for all lagged endogenous) -//datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, +//datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1, //istart_rmse=2, //start computing rmse from second observation (i.e. rmse does not inlude initial big error) //rmse=1, // do rmse analysis //); @@ -144,13 +144,13 @@ disp(' '); disp('Press ENTER to continue'); pause(5); // run this to generate posterior mode and Metropolis files if not yet done -estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2, +estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2, prefilter=1,mh_jscale=0.5,mh_replic=5000, mode_compute=4, mh_drop=0.6, nodisplay, bayesian_irf, filtered_vars, smoother) y_obs R_obs pie_obs dq de; // run this to produce posterior samples of filtered, smoothed and irf variables, if not yet done -//estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2,prefilter=1,mh_jscale=0.3, +//estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2,prefilter=1,mh_jscale=0.3, // mh_replic=0, mode_file=ls2003_mode, mode_compute=0, load_mh_file, bayesian_irf, // filtered_vars, smoother, mh_drop=0.6); @@ -163,7 +163,7 @@ disp('Press ENTER to continue'); pause(5); dynare_sensitivity(nodisplay, pprior=0,Nsam=2048,neighborhood_width=0.2, mode_file=ls2003_mode, // specifies the mode file where the mode and Hessian are stored -datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, +datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1, rmse=1); disp(' '); @@ -182,7 +182,7 @@ disp('RMSE ANALYSIS FOR MULTIVARIATE SAMPLE AT THE POSTERIOR MODE'); disp(' '); disp('Press ENTER to continue'); pause(5); dynare_sensitivity(nodisplay, mode_file=ls2003_mode, -datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, +datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1, pprior=0, stab=0, rmse=1, @@ -195,12 +195,12 @@ disp('THE LAST TWO CALLS COULD BE DONE TOGETHER'); disp('BY USING THE COMBINED CALL'); disp(' '); disp('dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file=ls2003_mode,') -disp('datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,') +disp('datafile=data_ca1.m,first_obs=8,nobs=79,prefilter=1,') disp('rmse=1, alpha2_rmse=0, alpha_rmse=0);') disp(' '); disp('Press ENTER to continue'); pause(5); //dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file=ls2003_mode, -//datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, +//datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1, //rmse=1, //alpha2_rmse=0, // no correlation analysis //alpha_rmse=0 // no Smirnov sensitivity analysis @@ -210,10 +210,10 @@ disp(' '); disp('RMSE ANALYSIS FOR POSTERIOR MCMC sample (ppost=1)'); disp('Needs a call to dynare_estimation to load all MH environment'); disp('Press ENTER to continue'); pause(5); -//estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2, mode_file=ls2003_mode, load_mh_file, +//estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2, mode_file=ls2003_mode, load_mh_file, // prefilter=1,mh_jscale=0.5,mh_replic=0, mode_compute=0, mh_drop=0.6); dynare_sensitivity(nodisplay, stab=0, // no need for stability analysis since the posterior sample is surely OK -datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, +datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1, rmse=1,ppost=1);