- Dev - dynare.org

Re: [DynareDev] [Trac] #9: Add derivatives of external functions in the model
by Trac Dynare 09 Jun '10

09 Jun '10

#9: Add derivatives of external functions in the model ---------------------------+------------------------------------------------ Reporter: sebastien | Owner: houtanb Type: feature | Status: accepted Priority: minor | Milestone: 4.2 Component: Preprocessor | Version: Resolution: | Keywords: ---------------------------+------------------------------------------------ Changes (by houtanb): * status: reopened => accepted -- Ticket URL: <https://www.dynare.org/trac/ticket/9#comment:9> Dynare <http://www.dynare.org> The Dynare project

1 0

[Trac] #100: Add exception for degenerate beta prior.
by Trac Dynare 09 Jun '10

09 Jun '10

#100: Add exception for degenerate beta prior. --------------------------+------------------------------------------------- Reporter: stepan | Owner: sebastien Type: bug | Status: new Priority: major | Milestone: 4.2 Component: Core M-files | Version: Keywords: | --------------------------+------------------------------------------------- If the declared prior expectation and standard deviation are equal to 0.5, the prior density is not defined because in this special case the parameters \alpha and \beta are zero. We should add an exception for this case. -- Ticket URL: <https://www.dynare.org/trac/ticket/100> Dynare <http://www.dynare.org> The Dynare project

1 3

Re: [DynareDev] Estimation DLL design - prior distributions
by Stéphane Adjemian 04 Jun '10

04 Jun '10

Hi George, The inverse gamma 2 density appears naturally in Bayesian statistics when considering conjugate priors in linear Gaussian models. The inverse gamma densities are quite common in Bayesian statistics (the references you are giving are not pertinent with respect to this literature). Furthermore the inverse gamma 2 is just a re-parametrization of the probability textbook (or common if you want) inverse gamma density. This is shown in the attached pdf. So I don't think we should provide the inverse gamma prior you are suggesting in Dynare. Best, Stéphane. On 2 June 2010 10:35, George Perendia <george(a)perendia.orangehome.co.uk>wrote: > Hi Michel > > Thanks for the pdf document . > > 2) My apologies - Having had a chance to see the full representation of > both inv gamma, 1 and 2 pdf in the document you sent I do not anymore > think there is a bug in the implementation - the implementation seem to > conform closely to the given representations. > > 1) I can not tell if the notes are incorrect but they do seem mutually > consistent, so that e.g. > > lpdfig1(x,s,n) = lpdfig2(x*x,s,n) = lpdfgam(1/(x*x),n/2,2/s) - > 2*log(x*x) = 0.5673 for my test > > On the other hand, both inv gamma 1 and 2 representations differ from what > I find the most common definition of the inverse gamma from various sources > > invGammaPDF(x,a,b) = b^a / gamma(a) * 1/x ^(a+1) * exp( -b/x) which, in > turn = GammaPDF(1/x,a,1/b) / x^2 > > (though, alternative representations of gammaPDF with inverse scale beta = > 1/b would give = GammaPDF(1/x,a,beta) / x^2 ) > > (e.g. see the enclosed document or internet reference such as > http://www.chemistrydaily.com/chemistry/Inverse-gamma_distribution > or even wikipedia > http://en.wikipedia.org/wiki/Inverse-gamma_distribution ) > > We may need to add another inv Gamma for that common representation. > > Best > George > > > On 01/06/2010 20:37, Michel Juillard wrote: > > Hi George, > > inverse gamma distributions are tricky because different authors use > different notations. I attach the notes that I had made in 2004. > Could you please tell me > 1) if you think that these notes are erroneous. In this case, please, > indicate the sources that you use. > 2) if the code currently implemented in Dynare Matlab conforms with these > notes. > > I'm leaving tomorrow for the Dynare Conference. I will be back this week > end and will try to look into the issue next week. > > All the best, > > Michel > > > On 06/01/2010 09:00 PM, George Perendia wrote: > > PS: On the other hand, using the commonly used inv gamma formula > > pdgInvGamma (x,n,s) = 1/(gamma(n) * s^n) * x ^(-n-1) exp (1/ (x*s) ) > > it appear that > > pdgInvGamma (x,n,s) = pdfGamma (1/x, n, s) / x^2 > > or, in logs, > > lpdfinvg (x,n,s) = -gammaln(n) - n*log(s) + (n+1)*log(1/x) - (1/x)/s > = [ -gammaln(n) - n*log(s) + (n-1)*log(1/x) - (1/x)/s ] - 2*log(x) = [ > lpdfgam(1/x,n,s) ] - 2*log(x) > > Best regards > > George > > On 01/06/2010 16:24, George Perendia wrote: > > Dear Michel > > thanks > > I have 2 issues: > > 1) According to few sources, > If x~InvGamma(a,b) , then 1/x ~Gamma(a,1/b) distribution > however, in our dynare inv gamma 2 (lpdfig2.m) we have different > definition: > X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) > (which is equivalent to > If x~InvGamma(a,b) , then 1/x ~Gamma(b/2,2/a) distribution) > > 2) I suspect there are a bugs in implementation for both invg2 and invg1 > even as per our specifications of invg2 aind invg1, i.e. > for invg2: we have > > ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) - > .5*(nu+2)*log(x(idx)) -.5*s./x(idx); > > and for e.g. s=2, n=5 and x= 0.5 as a result we have (incorrectly?): > > lpdfig2(x,s,n) = > -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n+2)*log(x) -.5*s./x > > ans = > > 0.1413 > > whilst > > lpdfgam(1/x,n/2,2/s) > > ans = > > -1.2450 > > and modified (corrected?) invgam2 with ...(nu - 2).which meets our spec.. > instead of (nu + 2)... > > K>> -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n-2)*log(x) -.5*s./x > > ans = > > -1.2450 > > thus, gives same result as lpdfgam(1/x,n/2,2/s). > > > Similarly, invg1 > > ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - > (nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ; > > lpdfig1(x,s,n)= 0.5673 = > lpdfig2(x*x,s,n)= 0.5673 > > whilst > lpdfgam(1/(x*x),n/2,2/s) > > ans = > > -2.2052 > > and modified invg1 without leading log(2) and with nu-2 instead nu+1 > > - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu-2)*log(x) - .5*s./(x.*x) > > ans = > > -2.2052 > > i.e. = lpdfgam(1/(x*x),n/2,2/s) > as per definition > > (... though I may have got something wrong) > > > If they are bugs, shall I log them and correct them? > > Best > george > > On 19/05/2010 20:29, Michel Juillard wrote: > > Dear George, > > thank you for the message. I agree that we should go with BOOST. > > 1) In the estimation DLL, we only need a uniform and a normal random > generator and the pdfs for the normal, beta, gamma, inverse gamma 1, inverse > gamma 2 and the uniform. You shouldn't try to implement more than that at > this stage. > > 2) for the inverse gamma 1 and 2 pdfs, you can directly translate the > current matlab function in ./matlab/distributions. We should add the source > of the algorithm if we can remember. > > 3) You need to look at the Matlab implementation, because we need the same > kind of generalizations: > - gamma is defined as greater than p3 rather than greater than 0 > - beta is defined over [p3; p4] instead of over [0,1] > > The relevant Matlab functions are set_prior.m and priordens.m > > All the best, > > Michel > > > > > On 05/19/2010 06:42 PM, George Perendia wrote: > > Hi Sebastien > > I started looking into priordens and pdfs as well as the relvant random > number generators (rng-s) for estimated parameters (*) and investigated two > options: C gsl and the C++ header BOOST as initially suggested by Stephane. > > Both provide for the most of the requirements but there are seeral > insufficiencies. > > gsl provides the necessary pair of pdf and random number generators > functions based on given uniform random number but not the uniform random > number generator - rng (or the pdf but that is trivial to add). > > On the other hand whilst BOOST provides a large variety of functions > including the pdfs and cdfs in a template form as well as a variety of > uniform and other distributions rngs and ways to combine them. Though it > does not seem to provide a rng for beta, the beta rng may be easily derived > using gamma rng and, e.g. the relation: > > "If *X* and *Y* are independently distributed Gamma<http://wiki/Gamma_distribution> > (*α*, *θ*) and Gamma(*β*, *θ*) respectively, then *X* / (*X* + *Y*) is > distributed Beta(*α*, *β*), " > > ...which is how gsl beta rng is implemented anyway: > > Neither gsl nor Boost however seem to provide pdfs or rngs for > inverse_gamma and that would need to be derived for either of them too. > > Whilst the gsl is simpler it is C whlst the BOOST is C++. Also, gsl would > require a library compilation before it can be used and it would also be an > additional package requirement for building Dynare (i.e. whilst we are > already using BOOST header library anyway). > > It seems that BOOST pdfs and rngs are better documented, provide a larger > variety of options for later extensions and fit better into our design > strategy whilst, as part of already used BOOST, they do not require any > extension on the overall maintenance and build requirements of Dynare. > > I wold be inclined therefore to suggest us to use BOOST but please do let > me know if you have any objections or other suggestions. > > If the performance is the issue, we can run a comparative test e.g. for > throwing 100000 beta rng-s and deriving their pdf-s.using gsl and BOOSt and > see the difference but that may a be bit excessive - just use of a different > uniform rng for generating specific distribution rng from a variety of those > available in BOOST can change performance dramatically. > > *) assuming the following distributions are required > beta_pdf | gamma_pdf | normal_pdf | uniform_pdf | inv_gamma_pdf | > inv_gamma1_pdf | inv_gamma2_pdf > > Best regards > George > > > On 18/05/2010 17:26, Sébastien Villemot wrote: > > Hi George, > > On Wed, May 12, 2010 at 09:36:26AM +0100, George Perendia wrote: > > > It was good to work together last week. > Aftert yesterday's ACES meeting it looks like I may not be engaged > full time on ACES for now, so, I could contiue on Estimation DLL, > 1-2 days a week for time being, e.g. on prior dens, optmiiser > integration, MCMC or something else. > Let me know what wold be best for you. > > > Sorry for my late answer. > > I am also happy with the progress that we made when you came to Paris. > > Unfortunately, since I have been working on other projects, I have not been able to finalize the likelihood DLL yet, but I am quite close to it. > > I am leaving for vacations until May 30, and I hope to finish the DLL by the first days of June. Note that I will still read and answer to emails during my vacations. > > Please feel free to advance on the implementation of other parts of the design. In particular, we need to decide how to implement the prior densities. We should rely on existing code rather than reimplementing all the probability density functions. My feeling is that using the GNU Scientific Library for that purpose is the easiest way, since all the PDFs are there and it is licensed under the GPL. Do you have any other option in mind ? > > Best, > > > > > > -- > > > > > > > > > > -- > > > > > > > -- > > > > > > > > > > -- > > Best regards > George > 07951415480 > > > > > -- Stéphane Adjemian CEPREMAP & Université du Maine Tel: (33)1-43-13-62-39

2 3

Re: [DynareDev] [Trac] #24: Add Markov-Switching
by Trac Dynare 01 Jun '10

01 Jun '10

#24: Add Markov-Switching ---------------------------+------------------------------------------------ Reporter: sebastien | Owner: Type: feature | Status: new Priority: major | Milestone: 4.2 Component: Core M-files | Version: Resolution: | Keywords: ---------------------------+------------------------------------------------ Comment(by sebastien): Note: the SWZ build system currently works on Linux and MacOS, but not Cygwin. The problem has to do with the GSL whichs creates conflicts with MATLAB libraries. -- Ticket URL: <https://www.dynare.org/trac/ticket/24#comment:1> Dynare <http://www.dynare.org> The Dynare project

1 0

[Trac] #104: Add option pruning
by Trac Dynare 31 May '10

31 May '10

#104: Add option pruning --------------------------+------------------------------------------------- Reporter: stepan | Owner: sebastien Type: bug | Status: new Priority: major | Milestone: 4.2 Component: Core M-files | Version: Keywords: | --------------------------+------------------------------------------------- Add option pruning in stoch_simul for second order simulations. The preprocessor should set options_.pruning equal to one if pruning option is used. The preprocessor should issue a warning if pruning and k_order_solver are used simultaneously (pruning is not available in dynare_simul_). -- Ticket URL: <https://www.dynare.org/trac/ticket/104> Dynare <http://www.dynare.org> The Dynare project

1 1

Re: [DynareDev] [Trac] #14: Create a database of old and new example MOD files
by Trac Dynare 30 Apr '10

30 Apr '10

#14: Create a database of old and new example MOD files ----------------------------+----------------------------------------------- Reporter: sebastien | Owner: Type: feature | Status: new Priority: major | Milestone: 4.2 Component: Documentation | Version: Resolution: | Keywords: ----------------------------+----------------------------------------------- Comment(by sebastien): Database started in commit 4a619f3db74b95aaca5c2a16a990d619df9975c4 -- Ticket URL: <https://www.dynare.org/trac/ticket/14#comment:3> Dynare <http://www.dynare.org> The Dynare project

1 0

[Trac] #46: Update examples on the website
by Trac Dynare 28 Apr '10

28 Apr '10

#46: Update examples on the website ---------------------------+------------------------------------------------ Reporter: sebastien | Owner: Type: bug | Status: new Priority: minor | Milestone: Component: Documentation | Version: Keywords: | ---------------------------+------------------------------------------------ Most MOD files examples distributed on the website don't work with Dynare 4. They were designed for Dynare 3. We regularly get posts on the forum by newcomers to Dynare who don't understand why these examples fail. We should update them to Dynare 4. -- Ticket URL: <https://www.dynare.org/trac/ticket/46> Dynare <http://www.dynare.org> The Dynare project

1 1

Estimation Parameters design
by George Perendia 28 Apr '10

28 Apr '10

Hi Sebastien (et al) On our EstimationModule design wiki page http://www.dynare.org/DynareWiki/EstimationModule we have few _id vectors indexing shocks and measurement errors correlation parameters into xparam1 such as.. * ncx_id ncx by 1 vector of integers.(for estimated structural shock correlations). . * ncn_id ncn by 1 vector of integers (estimated measurement error correlations). We have those *_id's indices into xparam1 but what indices are there into Q and H matrices respectively , including for pairs of (symmetric) coordinates on the lines of estim_params_.corrx(i,1); and k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1)); *used in DsgeLikelihood.m*? Shall we continue to use those as separate entities in estim_params_, including. estim_params_.var_exo and estim_params_.var_endo, or store separate new ones into the new EstimationModule params structure? -- Best regards George

2 4

New Dynare branches
by Sébastien Villemot 22 Apr '10

22 Apr '10

Hi everyone, I have recently created two new branches on the Dynare repository: * "gpm": a branch for keeping the current Dynare version used for the GPM (Global Projection Models) project * "dr1break": a topic branch to redesign dr1.m along the lines of: http://www.dynare.org/DynareWiki/BreakingDr1 This branch will be merged back in master when the job is finished. Unless you're explicitly involved in one of these two projects, you can safely ignore these branches. Best -- Sébastien Villemot

1 0

Kalman Filter accuracy
by George Perendia 21 Apr '10

21 Apr '10

Hi there I am implementing Kalman Filter in C++ There however appear that some apparently minor numerical discrepancies between Matlab and C++ versions of T and R matrices (i.e. the ghx and ghu matrices, and probably the other calculations too), are gradually creeping-up to a rather significan difference in the KF log-likelihood calculation and I am trying to figure out what can be done about it. The main impact seem to be on the inv(F) due to rather low values appearing in F. where F=ZPZ'+H resulting in a significant impact on the log likelihood and, as well, the next step P and system state forecast since both depend on the gain matrix K=P*Z' *invF and thus on invF; I.e. C++ F: F.data[0] = 0.00037206242498542022 F.data[1] =[2] = -0.00045917810077178716 F.data[3] = 0.0008672061055033099 and (rounded up) inverse of C++ F invF = [7756 4106.76; 4106.76 3327.6] (or using Matlab for inversion calculation invF = 1.0e+003 * 7.7560 4.1068 4.1068 3.3276) whilst Matlab Dynare F = 1.0e-003 * 0.3370 -0.4437 -0.4437 0.8383 and its iF = inv(F) = 1.0e+003 * 9.7907 5.1823 5.1823 3.9360 Then, e.g. Matlab likelihood factor v'*iF*v = 963.6491 whilst the C++ one is 771.88779, thus, discrepancies pretty much affecting the total LL calculation to say the least. It may also be a question which one is the more correct result. Any suggestions? PS: both, C++ and Matlab Dynare use the same C++<model>_dynamic.dll to calculate the initial jacobian but the jacobians differ. Best regards, George

2 2