From d38c4de498dacc6a188ca5c6fe476fdac6c9f93e Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Fri, 10 Dec 2010 21:34:09 +0100 Subject: [PATCH] estimation_dll: adding test for logposterior --- .../random_walk_metropolis_hastings_core.m | 285 ++++++ .../rawdata_euromodel_1.m | 967 ++++++++++++++++++ .../logposterior_dll_test/sweuromodel_dll.mod | 184 ++++ 3 files changed, 1436 insertions(+) create mode 100644 mex/sources/estimation/tests/logposterior_dll_test/random_walk_metropolis_hastings_core.m create mode 100644 mex/sources/estimation/tests/logposterior_dll_test/rawdata_euromodel_1.m create mode 100644 mex/sources/estimation/tests/logposterior_dll_test/sweuromodel_dll.mod diff --git a/mex/sources/estimation/tests/logposterior_dll_test/random_walk_metropolis_hastings_core.m b/mex/sources/estimation/tests/logposterior_dll_test/random_walk_metropolis_hastings_core.m new file mode 100644 index 000000000..afad53389 --- /dev/null +++ b/mex/sources/estimation/tests/logposterior_dll_test/random_walk_metropolis_hastings_core.m @@ -0,0 +1,285 @@ +function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab) +% PARALLEL CONTEXT +% This function contain the most computationally intensive portion of code in +% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for' +% cycle and are completely independent than suitable to be executed in parallel way. +% +% INPUTS +% o myimput [struc] The mandatory variables for local/remote +% parallel computing obtained from random_walk_metropolis_hastings.m +% function. +% o fblck and nblck [integer] The Metropolis-Hastings chains. +% o whoiam [integer] In concurrent programming a modality to refer to the differents thread running in parallel is needed. +% The integer whoaim is the integer that +% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the +% cluster. +% o ThisMatlab [integer] Allows us to distinguish between the +% 'main' matlab, the slave matlab worker, local matlab, remote matlab, +% ... Then it is the index number of this slave machine in the cluster. +% OUTPUTS +% o myoutput [struc] +% If executed without parallel is the original output of 'for b = +% fblck:nblck' otherwise a portion of it computed on a specific core or +% remote machine. In this case: +% record; +% irun; +% NewFile; +% OutputFileName +% +% ALGORITHM +% Portion of Metropolis-Hastings. +% +% SPECIAL REQUIREMENTS. +% None. + +% PARALLEL CONTEXT +% The most computationally intensive part of this function may be executed +% in parallel. The code sutable to be executed in parallel on multi core or cluster machine, +% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion. +% Then the DYNARE parallel package contain a set of pairs matlab functios that can be executed in +% parallel and called name_function.m and name_function_core.m. +% In addition in the parallel package we have second set of functions used +% to manage the parallel computation. +% +% This function was the first function to be parallelized, later other +% functions have been parallelized using the same methodology. +% Then the comments write here can be used for all the other pairs of +% parallel functions and also for management funtions. + + +% Copyright (C) 2006-2008,2010 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +if nargin<4, + whoiam=0; +end + + +global bayestopt_ estim_params_ options_ M_ oo_ + +% reshape 'myinputs' for local computation. +% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by: + +TargetFun=myinputs.TargetFun; +ProposalFun=myinputs.ProposalFun; +xparam1=myinputs.xparam1; +vv=myinputs.vv; +mh_bounds=myinputs.mh_bounds; +ix2=myinputs.ix2; +ilogpo2=myinputs.ilogpo2; +ModelName=myinputs.ModelName; +fline=myinputs.fline; +npar=myinputs.npar; +nruns=myinputs.nruns; +NewFile=myinputs.NewFile; +MAX_nruns=myinputs.MAX_nruns; +d=myinputs.d; +InitSizeArray=myinputs.InitSizeArray; +record=myinputs.record; +varargin=myinputs.varargin; + +% Necessary only for remote computing! +if whoiam + Parallel=myinputs.Parallel; + % initialize persistent variables in priordens() + priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ... + bayestopt_.p3,bayestopt_.p4,1); +end + +% (re)Set the penalty +bayestopt_.penalty = Inf; + +MhDirectoryName = CheckPath('metropolis'); + +options_.lik_algo = 1; +OpenOldFile = ones(nblck,1); +if strcmpi(ProposalFun,'rand_multivariate_normal') + n = npar; +elseif strcmpi(ProposalFun,'rand_multivariate_student') + n = options_.student_degrees_of_freedom; +end +% load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record'); +%%%% +%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains +%%%% + + +if any(isnan(bayestopt_.jscale)) + if exist([ModelName '_optimal_mh_scale_parameter.mat'])% This file is created by mode_compute=6. + load([ModelName '_optimal_mh_scale_parameter']) + proposal_covariance = d*Scale; + else + error('mh:: Something is wrong. I can''t figure out the value of the scale parameter.') + end +else + proposal_covariance = d*diag(bayestopt_.jscale); +end + + +jloop=0; + +for b = fblck:nblck, + jloop=jloop+1; + randn('state',record.Seeds(b).Normal); + rand('state',record.Seeds(b).Unifor); + if (options_.load_mh_file~=0) & (fline(b)>1) & OpenOldFile(b) + load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ... + '_blck' int2str(b) '.mat']) + x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)]; + logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)]; + OpenOldFile(b) = 0; + else + x2 = zeros(InitSizeArray(b),npar); + logpo2 = zeros(InitSizeArray(b),1); + end + if exist('OCTAVE_VERSION') || options_.console_mode + diary off + disp(' ') + elseif whoiam + % keyboard; + waitbarString = ['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']; + % waitbarTitle=['Metropolis-Hastings ',options_.parallel(ThisMatlab).ComputerName]; + if options_.parallel(ThisMatlab).Local, + waitbarTitle=['Local ']; + else + waitbarTitle=[options_.parallel(ThisMatlab).ComputerName]; + end + fMessageStatus(0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab)); + else, + hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']); + set(hh,'Name','Metropolis-Hastings'); + + end + isux = 0; + jsux = 0; + irun = fline(b); + j = 1; + while j <= nruns(b) + par = feval(ProposalFun, ix2(b,:), proposal_covariance, n); + if all( par(:) > mh_bounds(:,1) ) & all( par(:) < mh_bounds(:,2) ) + try + logpost = - feval(TargetFun, par(:),varargin{:}); + catch + logpost = -inf; + end + % testing logposterior DLL + [junk,logpost1] = logposterior(par(:),varargin{2},mexext); + if abs(logpost+logpost1) > 1e-10; + disp ([logpost -logpost1]) + end + else + logpost = -inf; + end + if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b)) + x2(irun,:) = par; + ix2(b,:) = par; + logpo2(irun) = logpost; + ilogpo2(b) = logpost; + isux = isux + 1; + jsux = jsux + 1; + else + x2(irun,:) = ix2(b,:); + logpo2(irun) = ilogpo2(b); + end + prtfrc = j/nruns(b); + if exist('OCTAVE_VERSION') || options_.console_mode + if mod(j, 10) == 0 + if exist('OCTAVE_VERSION') + printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j); + else + fprintf(' MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acceptance rate: %3.f \b%%\r', b, nblck, 100 * prtfrc, 100 * isux / j); + end + end + if mod(j,50)==0 & whoiam + % keyboard; + waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%%%', 100 * isux/j)]; + fMessageStatus(prtfrc,whoiam,waitbarString, '', options_.parallel(ThisMatlab)); + end + else + if mod(j, 3)==0 & ~whoiam + waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]); + elseif mod(j,50)==0 & whoiam, + % keyboard; + waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]; + fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab)); + end + end + + if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations + save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2'); + fidlog = fopen([MhDirectoryName '/metropolis.log'],'a'); + fprintf(fidlog,['\n']); + fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']); + fprintf(fidlog,' \n'); + fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']); + fprintf(fidlog,[' Acceptation rate......: ' num2str(jsux/length(logpo2)) '\n']); + fprintf(fidlog,[' Posterior mean........:\n']); + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']); + fprintf(fidlog,[' Minimum value.........:\n']);; + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']); + fprintf(fidlog,[' Maximum value.........:\n']); + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']); + fprintf(fidlog,' \n'); + fclose(fidlog); + jsux = 0; + if j == nruns(b) % I record the last draw... + record.LastParameters(b,:) = x2(end,:); + record.LastLogLiK(b) = logpo2(end); + end + % size of next file in chain b + InitSizeArray(b) = min(nruns(b)-j,MAX_nruns); + % initialization of next file if necessary + if InitSizeArray(b) + x2 = zeros(InitSizeArray(b),npar); + logpo2 = zeros(InitSizeArray(b),1); + NewFile(b) = NewFile(b) + 1; + irun = 0; + end + end + j=j+1; + irun = irun + 1; + end% End of the simulations for one mh-block. + record.AcceptationRates(b) = isux/j; + if exist('OCTAVE_VERSION') || options_.console_mode + if exist('OCTAVE_VERSION') + printf('\n'); + else + fprintf('\n'); + end + diary on; + elseif ~whoiam + close(hh); + end + record.Seeds(b).Normal = randn('state'); + record.Seeds(b).Unifor = rand('state'); + OutputFileName(jloop,:) = {[MhDirectoryName,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']}; +end% End of the loop over the mh-blocks. + + +myoutput.record = record; +myoutput.irun = irun; +myoutput.NewFile = NewFile; +myoutput.OutputFileName = OutputFileName; \ No newline at end of file diff --git a/mex/sources/estimation/tests/logposterior_dll_test/rawdata_euromodel_1.m b/mex/sources/estimation/tests/logposterior_dll_test/rawdata_euromodel_1.m new file mode 100644 index 000000000..cde44abc9 --- /dev/null +++ b/mex/sources/estimation/tests/logposterior_dll_test/rawdata_euromodel_1.m @@ -0,0 +1,967 @@ +C =[ + -7.4073 + -6.1860 + -6.5983 + -5.6088 + -5.0547 + -4.4774 + -3.8081 + -3.8425 + -2.4178 + -1.9835 + -1.0395 + -0.1583 + -0.0397 + 0.3505 + -0.1879 + -0.0067 + 0.0478 + -1.2247 + -1.4349 + -0.7973 + -0.0461 + 0.5844 + 1.1372 + 1.3801 + 1.8023 + 2.2972 + 2.0469 + 2.5435 + 2.8169 + 3.2007 + 2.6705 + 3.0518 + 3.2445 + 3.8443 + 3.8525 + 4.9494 + 4.2770 + 4.9532 + 5.1441 + 3.7124 + 3.9880 + 3.6926 + 2.6005 + 1.8679 + 1.9085 + 1.5563 + 1.2308 + 0.3264 + -0.2208 + -0.2483 + -0.4082 + -1.0315 + -1.6030 + -1.5499 + -1.3777 + -2.1675 + -2.5138 + -2.8820 + -2.6958 + -2.4719 + -1.9854 + -1.7954 + -2.2362 + -1.0595 + -0.8808 + -0.8548 + -1.2839 + -0.1363 + 0.2104 + 0.8810 + 0.3555 + 0.4766 + 1.3269 + 1.4506 + 1.4308 + 1.6263 + 1.9842 + 2.3948 + 2.8710 + 3.0177 + 2.9305 + 3.1739 + 3.7380 + 3.8285 + 3.3342 + 3.7447 + 3.7830 + 3.1039 + 2.8413 + 3.0338 + 0.3669 + 0.0847 + 0.0104 + 0.2115 + -0.6649 + -0.9625 + -0.7330 + -0.8664 + -1.4441 + -1.0179 + -1.2729 + -1.9539 + -1.4427 + -2.0371 + -1.9764 + -2.5654 + -2.8570 + -2.5842 + -3.0427 + -2.8312 + -2.3320 + -2.2768 + -2.1816 + -2.1043 + -1.8969 + -2.2388 + -2.1679 + -2.1172 +]; + +E =[ + 0.6263 + 0.7368 + 0.7477 + 1.0150 + 0.6934 + 0.4135 + 0.3845 + 0.2380 + 0.2853 + 0.5999 + 0.8622 + 1.2116 + 1.4921 + 1.5816 + 1.7259 + 1.6276 + 1.2422 + 0.8084 + 0.4710 + -0.3704 + -0.6427 + -0.5323 + -0.5562 + -0.3651 + -0.4356 + -0.7164 + -0.5816 + -0.4635 + -0.8456 + -0.9708 + -0.7138 + -0.7499 + -0.6941 + -0.6656 + -0.2912 + -0.1650 + 0.0774 + 0.2307 + 0.4484 + 0.4942 + 0.4653 + 0.2196 + 0.1736 + -0.1595 + -0.3918 + -0.4611 + -0.8493 + -0.7384 + -1.0604 + -1.2166 + -1.7187 + -1.6932 + -1.7830 + -1.7035 + -2.2079 + -2.3769 + -2.2511 + -2.1093 + -2.4638 + -2.4027 + -2.1313 + -1.9199 + -1.7941 + -1.4661 + -1.2269 + -1.0392 + -1.0725 + -0.7156 + -0.4778 + -0.4233 + -0.0409 + 0.1620 + 0.4280 + 0.5873 + 1.0323 + 1.3420 + 1.6902 + 2.0680 + 2.8219 + 3.2511 + 3.2930 + 3.5633 + 3.8992 + 3.6874 + 3.2849 + 3.1614 + 2.6221 + 2.5067 + 1.9223 + 1.1777 + 0.4483 + -0.0661 + -0.4424 + -0.9000 + -1.1478 + -1.2047 + -1.1412 + -1.2383 + -1.1048 + -0.9716 + -0.9287 + -1.0057 + -1.0827 + -1.0200 + -1.0072 + -1.1740 + -1.2809 + -1.1086 + -0.9866 + -0.8947 + -0.5875 + -0.2329 + 0.1493 + 0.4906 + 0.8400 + 1.0720 + 1.2648 + 1.5431 +]; + +I =[ + 2.6617 + 2.4325 + 1.9592 + 3.2530 + 2.9949 + 3.7918 + 4.7444 + 4.8289 + 5.5983 + 7.8923 + 9.4297 + 9.5010 + 10.0150 + 10.0413 + 9.6046 + 6.4766 + 5.9647 + 3.0114 + 0.5683 + -2.1226 + -2.1855 + -0.8329 + -1.5207 + -1.3419 + -1.7897 + -0.1476 + 0.4675 + -1.6516 + -1.5419 + -1.3050 + -1.2451 + -0.7815 + -0.7796 + -0.3612 + -2.4072 + 1.1162 + 1.1383 + 3.4132 + 5.0356 + 2.8016 + 2.1734 + 0.9366 + -0.7050 + -1.5021 + -2.9868 + -6.0237 + -6.2589 + -6.9138 + -8.2340 + -9.2589 + -9.2465 + -9.6988 + -9.7782 + -10.5645 + -10.7544 + -13.1583 + -12.2718 + -12.0131 + -13.5983 + -12.3579 + -10.9146 + -11.1572 + -12.4935 + -9.4393 + -8.5535 + -7.3723 + -10.0169 + -6.6088 + -5.2045 + -4.1024 + -2.8472 + -1.3139 + 0.0477 + 1.5629 + 3.6947 + 4.0327 + 4.1320 + 7.1400 + 9.1036 + 8.5609 + 7.6576 + 8.8022 + 8.9611 + 10.0871 + 9.4797 + 9.3964 + 10.0363 + 8.6340 + 6.6522 + 4.4471 + 0.2854 + -2.1879 + -2.9879 + -4.1021 + -2.7713 + -2.2281 + -1.2908 + -0.3250 + 0.6534 + 0.3942 + 0.3534 + -0.1532 + -1.7936 + 0.4909 + 0.3634 + 0.4290 + -0.9709 + 0.1942 + 0.6103 + 1.4426 + 2.7225 + 1.7525 + 3.2780 + 3.5985 + 4.9011 + 5.3312 + 6.4402 + 6.6529 +]; + +L =[ + 0.6263 + 0.7368 + 0.7477 + 1.0150 + 0.6934 + 0.4135 + 0.3845 + 0.2380 + 0.2853 + 0.5999 + 0.8622 + 1.2116 + 1.4921 + 1.5816 + 1.7259 + 1.6276 + 1.2422 + 0.8084 + 0.4710 + -0.3704 + -0.6427 + -0.5323 + -0.5562 + -0.3651 + -0.4356 + -0.7164 + -0.5816 + -0.4635 + -0.8456 + -0.9708 + -0.7138 + -0.7499 + -0.6941 + -0.6656 + -0.2912 + -0.1650 + 0.0774 + 0.2307 + 0.4484 + 0.4942 + 0.4653 + 0.2196 + 0.1736 + -0.1595 + -0.3918 + -0.4611 + -0.8493 + -0.7384 + -1.0604 + -1.2166 + -1.7187 + -1.6932 + -1.7830 + -1.7035 + -2.2079 + -2.3769 + -2.2511 + -2.1093 + -2.4638 + -2.4027 + -2.1313 + -1.9199 + -1.7941 + -1.4661 + -1.2269 + -1.0392 + -1.0725 + -0.7156 + -0.4778 + -0.4233 + -0.0409 + 0.1620 + 0.4280 + 0.5873 + 1.0323 + 1.3420 + 1.6902 + 2.0680 + 2.8219 + 3.2511 + 3.2930 + 3.5633 + 3.8992 + 3.6874 + 3.2849 + 3.1614 + 2.6221 + 2.5067 + 1.9223 + 1.1777 + 0.4483 + -0.0661 + -0.4424 + -0.9000 + -1.1478 + -1.2047 + -1.1412 + -1.2383 + -1.1048 + -0.9716 + -0.9287 + -1.0057 + -1.0827 + -1.0200 + -1.0072 + -1.1740 + -1.2809 + -1.1086 + -0.9866 + -0.8947 + -0.5875 + -0.2329 + 0.1493 + 0.4906 + 0.8400 + 1.0720 + 1.2648 + 1.5431 +]; + +PIE =[ + -1.0113 + -0.8305 + 0.2332 + -0.8746 + -0.7978 + -0.9220 + -0.2487 + -0.7749 + -0.5460 + -0.5347 + 0.5050 + -0.0334 + 0.6756 + 0.8791 + 0.7267 + 1.0997 + 1.1750 + 1.1927 + 0.4420 + 0.5357 + 0.0345 + 0.0196 + 0.3371 + 0.9379 + 1.2160 + 0.3393 + 0.5813 + 0.7410 + 0.3374 + 0.2616 + 0.4025 + 0.4799 + 0.5981 + -0.1523 + 0.4458 + 0.2182 + 0.9793 + 0.7562 + 1.0064 + 0.8203 + 0.6966 + 0.3352 + 0.6581 + 0.6111 + 0.9833 + 1.1991 + 0.9562 + 0.3868 + 0.2939 + 0.2471 + 0.8331 + 0.0715 + 0.3910 + 0.3301 + 0.2547 + -0.2702 + -0.2998 + -0.1953 + -0.2293 + -0.3284 + 0.0480 + -0.0374 + 0.3253 + -0.3434 + -0.3892 + -0.7178 + -0.4758 + -0.6794 + -0.8505 + -0.3512 + -0.4436 + -0.5101 + -0.4574 + -0.2696 + -0.1047 + -0.5745 + -0.2989 + -0.0063 + 0.0088 + -0.1184 + -0.1506 + -0.4073 + 0.2674 + 0.2896 + 0.0669 + 0.1166 + -0.1699 + -0.2518 + -0.0562 + -0.3269 + -0.0703 + -0.1046 + -0.4888 + -0.3524 + -0.2485 + -0.5870 + -0.4546 + -0.3970 + -0.2353 + -0.0352 + -0.2171 + -0.3754 + -0.4322 + -0.4572 + -0.4903 + -0.4518 + -0.6435 + -0.6304 + -0.4148 + -0.2892 + -0.4318 + -0.6010 + -0.4148 + -0.4315 + -0.3531 + -0.8053 + -0.4680 + -0.4263 +]; + +R =[ + -1.0750 + -1.1540 + -1.3682 + -1.4569 + -1.3490 + -1.4011 + -1.6486 + -1.6968 + -1.6976 + -1.2567 + -1.1392 + -0.7783 + -0.3021 + -0.0435 + 0.0066 + -0.0043 + 0.1029 + -0.0628 + -0.5358 + -0.9627 + -1.1079 + -1.0918 + -0.9966 + -0.6223 + -0.3616 + -0.2711 + -0.0997 + -0.2810 + -0.3710 + -0.3167 + -0.5301 + -0.5826 + -0.3194 + -0.2713 + -0.5287 + -0.2432 + 0.1098 + 0.5349 + 0.7094 + 0.8415 + 0.6226 + 0.7376 + 0.9316 + 1.4370 + 1.5853 + 1.4267 + 1.1783 + 1.2046 + 0.9689 + 0.7918 + 0.6315 + 0.5950 + 0.6853 + 0.7171 + 0.5887 + 0.4873 + 0.4027 + 0.3489 + 0.2934 + 0.3060 + 0.1741 + 0.0348 + 0.0771 + -0.1005 + -0.1518 + -0.1104 + -0.0681 + -0.0059 + 0.0256 + 0.0404 + -0.1721 + -0.2002 + 0.0015 + 0.1249 + 0.3738 + 0.4320 + 0.5579 + 0.8186 + 0.8727 + 0.7356 + 0.7243 + 0.8635 + 0.9058 + 0.7656 + 0.7936 + 0.8631 + 0.9074 + 0.9547 + 1.2045 + 1.0850 + 0.9178 + 0.5242 + 0.3178 + 0.1472 + 0.0227 + -0.0799 + -0.0611 + -0.0140 + 0.1132 + 0.1774 + 0.0782 + 0.0436 + -0.1596 + -0.2691 + -0.2895 + -0.3791 + -0.4020 + -0.4166 + -0.4037 + -0.3636 + -0.4075 + -0.4311 + -0.4470 + -0.5111 + -0.6274 + -0.7261 + -0.6974 + -0.5012 +]; + +W =[ + -14.8791 + -13.2300 + -13.5037 + -13.0249 + -11.2546 + -10.0148 + -8.8586 + -8.5739 + -7.7851 + -6.7136 + -5.5878 + -4.6881 + -3.8039 + -3.0366 + -2.7342 + -1.3135 + -0.7387 + -0.1131 + -0.2769 + 0.8696 + 1.8855 + 2.3667 + 2.4942 + 3.2049 + 3.9682 + 5.1500 + 4.7047 + 4.7827 + 5.3377 + 5.6614 + 5.2813 + 5.2967 + 5.5175 + 6.1526 + 5.6627 + 6.0694 + 6.5824 + 6.9032 + 6.7849 + 6.6896 + 6.6201 + 6.9933 + 5.8959 + 6.7419 + 6.9999 + 6.4009 + 5.5083 + 5.1054 + 5.2813 + 4.5790 + 3.9589 + 3.8599 + 3.8978 + 2.7957 + 3.2480 + 1.4634 + 1.9219 + 1.8398 + 1.9279 + 1.8316 + 1.6092 + 1.2741 + 0.2031 + -0.0236 + -0.1004 + -0.3034 + -1.0273 + -0.2205 + 0.0458 + 0.2386 + -0.0977 + -0.3145 + -0.1416 + -0.7009 + -0.9082 + -0.8802 + -0.5644 + -0.5852 + -0.5346 + 0.0652 + 0.1301 + 0.3444 + -0.3592 + 0.8096 + 0.9644 + 1.0289 + 1.2781 + 1.2298 + 2.2134 + 2.0808 + 0.4925 + 0.6506 + 0.5531 + 0.2456 + -0.5351 + -0.8183 + -0.8967 + -0.7268 + -1.0738 + -1.2844 + -1.4338 + -1.6995 + -1.7085 + -2.2889 + -2.1018 + -2.4273 + -2.4609 + -2.1407 + -2.3847 + -3.1689 + -4.5581 + -4.1027 + -4.2436 + -4.8836 + -5.9660 + -4.9971 + -5.2386 + -5.6618 +]; + +Y =[ + -4.9347 + -4.6205 + -5.2198 + -4.5937 + -3.8015 + -3.6643 + -2.7239 + -2.7524 + -2.0634 + -1.0112 + 0.0530 + 0.7623 + 1.7927 + 2.1486 + 2.4866 + 2.1456 + 2.1671 + -0.0254 + -1.6716 + -1.9673 + -1.6109 + -1.0292 + -0.1222 + 0.7329 + 1.1234 + 2.0603 + 1.7998 + 1.4820 + 1.1732 + 1.6424 + 1.5382 + 2.1399 + 2.0127 + 2.7210 + 2.4966 + 3.5249 + 3.6237 + 4.2011 + 4.5634 + 3.3442 + 2.7761 + 1.9812 + 1.3779 + 1.4616 + 1.3029 + 0.7594 + 0.3695 + 0.0832 + -0.8118 + -1.4557 + -1.4850 + -1.2346 + -1.5696 + -1.3785 + -0.7682 + -2.0308 + -1.7778 + -1.7801 + -2.1711 + -1.7469 + -1.3413 + -1.3352 + -2.4390 + -1.2125 + -1.1695 + -1.0891 + -2.4753 + -1.3503 + -0.9412 + -0.1470 + 0.0026 + 0.1108 + 0.6890 + 1.3520 + 1.6018 + 2.0667 + 1.7625 + 2.6658 + 3.4048 + 3.2507 + 3.4251 + 3.2174 + 3.1903 + 3.3396 + 3.1358 + 2.8625 + 3.3546 + 2.4609 + 1.9534 + 0.9962 + -0.7904 + -1.1672 + -1.2586 + -1.3593 + -1.3443 + -0.9413 + -0.6023 + -0.4516 + -0.5129 + -0.8741 + -1.0784 + -1.4091 + -1.3627 + -1.5731 + -1.6037 + -1.8814 + -2.1482 + -1.3597 + -1.1855 + -1.1122 + -0.8424 + -0.9747 + -1.1385 + -1.4548 + -1.4284 + -1.4633 + -1.0621 + -0.7871 +]; diff --git a/mex/sources/estimation/tests/logposterior_dll_test/sweuromodel_dll.mod b/mex/sources/estimation/tests/logposterior_dll_test/sweuromodel_dll.mod new file mode 100644 index 000000000..999485ee8 --- /dev/null +++ b/mex/sources/estimation/tests/logposterior_dll_test/sweuromodel_dll.mod @@ -0,0 +1,184 @@ +//options_.usePartInfo=1; + +var MC E EF R_KF QF CF IF YF LF PIEF WF RF R_K Q C I Y L PIE W R EE_A PIE_BAR EE_B EE_G EE_L EE_I KF K one BIGTHETA; + +varexo E_A E_B E_G E_L E_I ETA_R E_PIE_BAR ETA_Q ETA_P ETA_W ; + +parameters +xi_e lambda_w alpha czcap beta phi_i tau sig_c hab ccs cinvs phi_y gamma_w xi_w gamma_p xi_p sig_l r_dpi +r_pie r_dy r_y rho rho_a rho_pb rho_b rho_g rho_l rho_i LMP ; + + + +alpha=.30; +beta=.99; +tau=0.025; +ccs=0.6; +cinvs=.22; //% alpha*(tau+ctrend)/R_K R_K=ctrend/beta-1+tau +lambda_w = 0.5; +phi_i= 6.771; +sig_c= 1.353; +hab= 0.573; +xi_w= 0.737; +sig_l= 2.400; +xi_p= 0.908; +xi_e= 0.599; +gamma_w= 0.763; +gamma_p= 0.469; +czcap= 0.169; +phi_y= 1.408; +r_pie= 1.684; +r_dpi= 0.14; +rho= 0.961; +r_y= 0.099; +r_dy= 0.159; +rho_a= 0.823; +rho_b= 0.855; +rho_g= 0.949; +rho_l= 0.889; +rho_i= 0.927; +rho_pb= 0.924; +LMP = 0.0 ; //NEW. + +model(linear, use_dll); + CF = (1/(1+hab))*(CF(1)+hab*CF(-1))-((1-hab)/((1+hab)*sig_c))*(RF-PIEF(1)-EE_B) ; + 0 = alpha*R_KF+(1-alpha)*WF -EE_A ; + PIEF = 0*one; + IF = (1/(1+beta))* (( IF(-1) + beta*(IF(1)))+(1/phi_i)*QF)+0*ETA_Q+EE_I ; + QF = -(RF-PIEF(1))+(1-beta*(1-tau))*((0+czcap)/czcap)*R_KF(1)+beta*(1-tau)*QF(1) +0*EE_I ; + KF = (1-tau)*KF(-1)+tau*IF(-1) ; + YF = (ccs*CF+cinvs*IF)+EE_G ; + + YF = 1*phi_y*( alpha*KF+alpha*(1/czcap)*R_KF+(1-alpha)*LF+EE_A ) ; + WF = (sig_c/(1-hab))*(CF-hab*CF(-1)) + sig_l*LF - EE_L ; + LF = R_KF*((1+czcap)/czcap)-WF+KF ; + EF = EF(-1)+EF(1)-EF+(LF-EF)*((1-xi_e)*(1-xi_e*beta)/(xi_e)); + + C = (hab/(1+hab))*C(-1)+(1/(1+hab))*C(1)-((1-hab)/((1+hab)*sig_c))*(R-PIE(1)-EE_B) ; + I = (1/(1+beta))* (( I(-1) + beta*(I(1)))+(1/phi_i)*Q )+1*ETA_Q+1*EE_I ; + Q = -(R-PIE(1))+(1-beta*(1-tau))*((0+czcap)/czcap)*R_K(1)+beta*(1-tau)*Q(1) +EE_I*0+0*ETA_Q ; + K = (1-tau)*K(-1)+tau*I(-1) ; + Y = (ccs*C+cinvs*I)+ EE_G ; + Y = phi_y*( alpha*K+alpha*(1/czcap)*R_K+(1-alpha)*L ) +phi_y*EE_A ; + PIE = (1/(1+beta*gamma_p))* + ( + (beta)*(PIE(1)) +(gamma_p)*(PIE(-1)) + +((1-xi_p)*(1-beta*xi_p)/(xi_p))*(MC) + ) + ETA_P ; + + MC = alpha*R_K+(1-alpha)*W -EE_A; + W = (1/(1+beta))*(beta*W(+1)+W(-1)) + +(beta/(1+beta))*(PIE(+1)) + -((1+beta*gamma_w)/(1+beta))*(PIE) + +(gamma_w/(1+beta))*(PIE(-1)) + -(1/(1+beta))*(((1-beta*xi_w)*(1-xi_w))/(((1+(((1+lambda_w)*sig_l)/(lambda_w))))*xi_w))*(W-sig_l*L-(sig_c/(1-hab))*(C-hab*C(-1))+EE_L) + +ETA_W; + L = R_K*((1+czcap)/czcap)-W+K ; + +// R = r_dpi*(PIE-PIE(-1)) +// +(1-rho)*(r_pie*(PIE(-1)-PIE_BAR)+r_y*(Y-YF)) +// +r_dy*(Y-YF-(Y(-1)-YF(-1))) +// +rho*(R(-1)-PIE_BAR) +// +PIE_BAR +// +ETA_R; + + + R = + +r_dpi*(PIE-PIE(-1)) + + +(1-rho)*(r_pie*(BIGTHETA)+r_y*(Y-YF)) + +r_dy*(Y-YF-(Y(-1)-YF(-1))) + +rho*(R(-1)-PIE_BAR) + +PIE_BAR + +ETA_R; + + + E = E(-1)+E(1)-E+(L-E)*((1-xi_e)*(1-xi_e*beta)/(xi_e)); + + + EE_A = (rho_a)*EE_A(-1) + E_A; + PIE_BAR = rho_pb*PIE_BAR(-1)+ E_PIE_BAR ; + EE_B = rho_b*EE_B(-1) + E_B ; + EE_G = rho_g*EE_G(-1) + E_G ; + EE_L = rho_l*EE_L(-1) + E_L ; + EE_I = rho_i*EE_I(-1) + E_I ; + one = 0*one(-1) ; + + LMP*BIGTHETA(1) = BIGTHETA - (PIE(-1) - PIE_BAR) ; + +end; + + +shocks; +var E_A; stderr 0.598; +var E_B; stderr 0.336; +var E_G; stderr 0.325; +var E_I; stderr 0.085; +var E_L; stderr 3.520; +var ETA_P; stderr 0.160; +var ETA_W; stderr 0.289; +var ETA_R; stderr 0.081; +var ETA_Q; stderr 0.604; +var E_PIE_BAR; stderr 0.017; +end; + +//stoch_simul(irf=20) Y C PIE R W R_K L Q I K ; + +// stoch_simul generates what kind of standard errors for the shocks ? + +//steady; +//check; +//stoch_simul(periods=200,irf=20,simul_seed=3) Y C PIE MC R W R_K E L I ; + +//datatomfile('ddd',[]); + +// new syntax + +estimated_params; +// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE +// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF +stderr E_A,0.543,0.01,4,INV_GAMMA_PDF,0.4,2; +stderr E_PIE_BAR,0.072,0.001,4,INV_GAMMA_PDF,0.02,10; +stderr E_B,0.2694,0.01,4,INV_GAMMA_PDF,0.2,2; +stderr E_G,0.3052,0.01,4,INV_GAMMA_PDF,0.3,2; +stderr E_L,1.4575,0.1,6,INV_GAMMA_PDF,1,2; +stderr E_I,0.1318,0.01,4,INV_GAMMA_PDF,0.1,2; +stderr ETA_R,0.1363,0.01,4,INV_GAMMA_PDF,0.1,2; +stderr ETA_Q,0.4842,0.01,4,INV_GAMMA_PDF,0.4,2; +stderr ETA_P,0.1731,0.01,4,INV_GAMMA_PDF,0.15,2; +stderr ETA_W,0.2462,0.1,4,INV_GAMMA_PDF,0.25,2; +rho_a,.9722,.1,.9999,BETA_PDF,0.85,0.1; +rho_pb,.85,.1,.999,BETA_PDF,0.85,0.1; +rho_b,.7647,.1,.99,BETA_PDF,0.85,0.1; +rho_g,.9502,.1,.9999,BETA_PDF,0.85,0.1; +rho_l,.9542,.1,.9999,BETA_PDF,0.85,0.1; +rho_i,.6705,.1,.99,BETA_PDF,0.85,0.1; +phi_i,5.2083,1,15,NORMAL_PDF,4,1.5; +sig_c,0.9817,0.25,3,NORMAL_PDF,1,0.375; +hab,0.5612,0.3,0.95,BETA_PDF,0.7,0.1; +xi_w,0.7661,0.3,0.9,BETA_PDF,0.75,0.05; +sig_l,1.7526,0.5,5,NORMAL_PDF,2,0.75; +xi_p,0.8684,0.3,0.95,BETA_PDF,0.75,0.05; +xi_e,0.5724,0.1,0.95,BETA_PDF,0.5,0.15; +gamma_w,0.6202,0.1,0.99,BETA_PDF,0.75,0.15; +gamma_p,0.6638,0.1,0.99,BETA_PDF,0.75,0.15; +czcap,0.2516,0.01,2,NORMAL_PDF,0.2,0.075; +phi_y,1.3011,1.001,2,NORMAL_PDF,1.45,0.125; +r_pie,1.4616,1.2,2,NORMAL_PDF,1.7,0.1; +r_dpi,0.1144,0.01,0.5,NORMAL_PDF,0.3,0.1; +rho,0.8865,0.5,0.99,BETA_PDF,0.8,0.10; +r_y,0.0571,0.01,0.2,NORMAL_PDF,0.125,0.05; +r_dy,0.2228,0.05,0.5,NORMAL_PDF,0.0625,0.05; +end; + +varobs Y C I E PIE W R; + +//estimation(datafile=rawdata_euromodel_1,presample=40, first_obs=1, nobs=118, lik_init=2, mode_compute=1,mh_replic=0); +estimation(datafile=rawdata_euromodel_1,presample=40, first_obs=1, nobs=118,mh_jscale=0.2,mh_replic=1000 +//,mode_compute=0,mode_file=sweuromodel_dll_mode +); + + +//stoch_simul(periods=200,irf=20,simul_seed=3) Y C PIE R W R_K L Q I K ; + -- GitLab