June 2009 - Dev - dynare.org

Re: [DynareDev] Quasi triangular matrices in Kalman filters
by Michel Juillard 24 Jun '09

24 Jun '09

Andrea and George, 1) please write a standalone test/timing of the QT code so that we can profile it using standard tools 2) Compare with call to dgmev() using Lapack + optimized Blas (possibly from a Matlab distribution). 3) Upload the code to SVN, so we can test it on other machines. The huge variability between 2 runs reported by George, may be due to Windows and usually is less important under Linux. All the best, Michel G. Perendia wrote: > Dear Andrea > > Thanks for the new libraries. > > I run some initial performance tests today for the simple T*a matrix*vector > multiplication with 2 different QT matrices sizes but in summary this is > what I am, getting: > > 10000 iteration loop with 100x100 random QT matrix (from qz decomposition) > and a vector: > > 1st & 2nd run (after restart) > > Native matlab matrix multiplication in a loop > > Ta_time = 0.3010 & 0.6610 > > Calling dgemv() using Sylv Vector and General Matrix is faster than Matlab > loop: > > GMcppTaInnrLoop_time = 0.1600 & 0.3310 > > Calling QT f90 library using of Sylv Vector and General Matrix: > > QTcppTaInnerLoop_time = 8.5730 & 20.5300 > > Calling QT f90 library without use of Sylv Vector and General Matrix but > using only pure C/C++ double arrays is only marginally faster: > > QTcpp_noSylv_TaInnerLoop_time = 8.4420 > > > > 1000 loop with 10x10 random QT matrix and vector: > > For a 10x10 matrix, calling QT f90 library takes about twice the time Matlab > loop does but dgemv is still faster. > > Matlab: 0.0400 > > GMcppTaInnrLoop_time = 0.0300 > > QTcppTaInnerLoop_time = 0.0800 > > > > It is, however, possible that the MinGW f95 I am using is not the best > optimising compiler that can be used and/or that tests for PTP', which I am > planning to do next, may be better.. > > > > What are your thoughts? Do you think that we may be able to improve > performance of this multiplication somehow. > > I wander if making many cross-language calls may be rather detrimental and > that we may improve performance if we reduce this high level of > modularisation and calling, e.g. by using a higher level subroutine that > will perform all operations within f90, passing back only the final Ta? > > > > NOTES: > > After a restart, Matlab appears to be much slower than later! > > Also, matlab multiplication reports both, the real and the imaginary part of > the result which appear complex but the real part matches QT and dgemv > outputs.. > > > > Best regards > > > George > artilogica(a)btconnect.com > > ----- Original Message ----- > From: "Michel Juillard" <michel.juillard(a)ens.fr> > To: "andrea pagano" <pagano.andrea(a)gmail.com> > Cc: "G. Perendia" <george(a)perendia.orangehome.co.uk> > Sent: Friday, June 19, 2009 1:33 PM > Subject: Re: Quasi triangular matrices in Kalman filters > > > >> Thanks Andrea >> >> amities >> >> Michel >> >> andrea pagano wrote: >> >>> Hi all >>> I would go for subroutines. I will do it over the weekend while >>> looking at other possibilities fortran pointers. >>> >>> >>> Best >>> >>> Andrea >>> >>> >>> On Fri, Jun 19, 2009 at 10:04 AM, G. >>> Perendia<george(a)perendia.orangehome.co.uk> wrote: >>> >>> >>>> Dear Andrea >>>> >>>> Problem: >>>> >>>> I have encountered a problem integrating KalmanFilter with the >>>> f90 QT library - passing the QT result arrays back to C++. >>>> >>>> QT Fortran routines have been written in standard Fortran FUNCTION >>>> > format, > >>>> (i.e., not SUBROUTINE), so that they are returning double or single >>>> dimensional array (they are named by), by value ( not reference). >>>> > However, > >>>> as it appears, only simple, single variables seems can be passed from >>>> Fortran >>>> FUNCTIONs back to C++ (e.g. INT or REAL). >>>> >>>> On the other hand, NAG, BLAS and LAPACK routines have all been written >>>> > as > >>>> Fortran SUBROUTINEs and they can be integrated with C more easily - >>>> > they > >>>> receive parameters >>>> and return their results through the variables passed as calling >>>> > parameters, > >>>> by references. >>>> >>>> For example, dgemv.f from BLAS library gets Y by reference and returns >>>> modified Y passed as calling parameter reference. >>>> >>>> SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) >>>> .... >>>> * Y - DOUBLE PRECISION array of DIMENSION at least >>>> ... >>>> * Before entry .... the incremented array Y >>>> * must contain the vector y. On exit, Y is overwritten by the >>>> * updated vector y. >>>> .... >>>> >>>> Poss. Solutions: >>>> >>>> I could not find any references on how to get arrays from Fortran >>>> > FUCTION as > >>>> return value back to C - do you or anyone around you know how to do it, >>>> > if > >>>> at all possible? In any case, passing array by value is also not >>>> > recommended > >>>> as rather un-economical, especially for larger matrices. >>>> >>>> One alternative way I can think of is less explored option of returning >>>> Fortran pointer to the resulting array from the QT functions instead >>>> > of the > >>>> array by value >>>> and I think I can work one of that out but suggestions are more than >>>> welcome. >>>> >>>> I can see few options: >>>> >>>> a) to rewrite QT library as SUBROUTINE instead FUCTION routines, or, >>>> >>>> b) try to use Fortran pointers and, if we can >>>> then also rewrite QT library to return pointers, or, >>>> >>>> c) write 3 or more high level f90 shell SUBROUTINES calling the >>>> > existing and > >>>> unmodified QT functions and performing the all operations needed to >>>> construct the resulting Ta and >>>> TPT' (for both cases 1 and 2) instead of doing low-level QT >>>> > manipulation in > >>>> C++. >>>> >>>> This way QT library need not be changed and those new SUBROUTINEs will >>>> > also > >>>> act as interface with C++. I think this is a more productive and >>>> > optimal > >>>> alternative of the three since those combination utilities would have >>>> > to be > >>>> written anyway, except it seem to be easier to do that now in Fortran >>>> > than > >>>> in C++. >>>> >>>> If you like and/or are busy, I can by Monday develop the Ta and the >>>> > first > >>>> case of TPT' SUBROUTINES whilst the second case may need more thinking >>>> > and > >>>> more granular approach to take advantage of multiple processors. >>>> >>>> >>>> Please let me know your thoughts on this issue and, whether if you have >>>> > time > >>>> to make the needed changes or additions in the f90 files. >>>> >>>> Best regards >>>> >>>> George >>>> artilogica(a)btconnect.com >>>> >>>> ----- Original Message ----- >>>> From: "andrea pagano" <pagano.andrea(a)gmail.com> >>>> To: <george(a)perendia.orangehome.co.uk> >>>> Cc: "Michel Juillard" <Michel.Juillard(a)ens.fr> >>>> Sent: Monday, June 01, 2009 7:47 PM >>>> Subject: Quasi triangular matrices in Kalman filters >>>> >>>> >>>> >>>> >>>>> Hi all >>>>> >>>>> I am sending you a set of Fortran routines to calculate the matricial >>>>> expression in Kalman filter together with some explenations. >>>>> >>>>> Hope they can be a starting point in optimizing the overall >>>>> > computation > >>>>> Best >>>>> >>>>> Andrea >>>>> >>>>> >>>>> -- >>>>> Andrea Pagano >>>>> via Veratti VARESE >>>>> tel. +3903321691261 >>>>> cell.+393403804397 >>>>> >>>>> >>>>> >>>> >>> >>> >>> >> > > > >

2 1

Re: [DynareDev] hessian from csminwel
by Michel Juillard 24 Jun '09

24 Jun '09

Dear George, optimziation routines such as csminwel (or fmincon) need only a rough Hessian, and use brute force to keep it definite positive. It is not appropriate to compute the standard error of the posterior mode, that is why we recompute it using a 7-point formula, but, of course, you need to really be at the correct mode to obtain a positive definite hessian. I imagine that in some cases the hessian returned by the optimization routine can do as variance matrix for the proposal, but in other cases, it can't. In doubt, we force the user the refine his/her mode before going into Metropolis. Best, Michel G. Perendia wrote: > Dear Michel > > Working on Part Info integration we are experiencing frequent problem > with models that run only to end of the first, the ML /mode stage > and fail to run MCMC because they fail to calculate chol(vv), > this due to the hessian of the estimated parameters not being positive > definite - this problem is hampering PI tests: > > ??? Error using ==> chol > Matrix must be positive definite. > Error in ==> metropolis_hastings_initialization at 52 > d = chol(vv); > > Why is Dynare using the additional parameter hessian calculation (i.e > using hessian()) in dynare_estimation_1.m when the last one (already > inverted) returned from csminwel (the default gradient optimiser, or > from fmincon) and passes chol(), may be used for MCMC stage > instead (i.e. as is the situation for the optimiser cases 5 - newrat & > 6)? > > What is disadvantage of using hessian_csminwel instead? (I > do understand though that some of the other optimisers do not produce > the Hessian and such calculation is then needed) > > Best regards > > George > artilogica(a)btconnect.com <mailto:artilogica@btconnect.com>

1 0

Release of Dynare 4.0.4
by Sébastien Villemot 20 Jun '09

20 Jun '09

Hi all, I made the packages for Dynare 4.0.4. The announce is as follows: Changes: * Fixed forecasting after stoch_simul * Fixed crash of the smoother when only one exogenous variable * Fixed issue with compilation of MEX files under 64 bits platforms with Matlab R2009a * Fixed name clash when MOD file is called "ramsey.mod" and computes optimal policy * Fixed USE_DLL option at order 2 * Updates to the reference manual in sections on model resolution and estimation Release notes: * for users of Octave under Windows, you must upgrade your Octave installation to 3.0.5 release, for compatibility with the new packages (see http://www.dynare.org/DynareWiki/DynareOctave ) * Linux packages are available for Ubuntu "Jaunty Jackalope" (9.04), Ubuntu "Intrepid Ibex" (8.10), Ubuntu "Hardy Heron" (8.04) and Debian "Lenny" (5.0) as usual (see http://www.dynare.org/DynareWiki/InstallOnDebianOrUbuntu ) @Michel: the stuff you need for the web site is in /home/villemot on Pythie. You can find there: * the Windows package * the source in .tar.gz * the reference manual in PDF and HTML (more recent than those currently on the web site) Best, -- Sébastien Villemot

1 0

Re: [DynareDev] Kalman Filter
by Michel Juillard 18 Jun '09

18 Jun '09

Thanks George, very interesting. Please proceed with simplifying the KF code. All the best, Michel G. Perendia wrote: > Dear Michel > > I already have a stand-alone exe test I used last week (and uploaded it > today too) and I run it through gprof earlier today (though, after wasting > some time last week trying to use a reportedly sophisticated profiler - > CodeAnalyst from AMD - which, so far, I could not make to work at all). > > The other profiling tool I initially used (and reported) last week ("Very > Sleepy") is an "external", not getting into the code details itself but > could be run attached externally to either the stand-alone exe test or the > Matlab running DLL thread too, and that reported for both spending a lot > (~40%) time in ztrsv solver but, at one snapshot, also pointed to a lot of > time spent in dtrsv and GeneralMatrix copy - 50% each. On the contrast, > gprof is more internal, higher resolution profiler and it puts the load > weight (>10%) on housekeeping functions but does not even mention calls to > external library BLAS functions such as dtrsv and ztrsv. > > The both profiling tools, however, seem to confirm what my early code > inspection concluded too: a very high, use of (not very productive) > General Matrix copy constructor -(e.g. the C++ kalman filter stores copy of > few of the main time variant system matrices: F, P and an intermediate L > for each step of the time series evaluation and also creates a copy of the > input T, H and Z at each step as if they may be time invariant too although > this would not be the case for a simple, no-diffuse KF without missing > observations). > > This then resulted in high % of time in GeneralMatrix copy() function, > (which is called by the copy-constructor explicitly) - i.e. as reported by > the both profiling programs: Sleepy: up to 50% at one snapshot point whilst > gprof gives it the 1st rank with 11% on its own or 27.3% of total with its > children. > > The copy() function is followed by the utility functions such as two > varieties (const and non-const) of Vector indexing [] operator and by the > const and variable varieties of GeneralMatrix::get() elements that utilise > the previous Vector indexing[] operator and are themselves directly called > from the heavily used GenaralMatrix copy function among the rest. > > According to gprof, the above high burden of copy constructor and the > related functions are only then followed by the productive functions such as > PLUFact::multInvRight matrix multiplication with inversion (used for > inversion of the F matrix), the GeneralMatrix constructors and the > GeneralMatrix::gemm() - a general matrix multiplication(itself calling > BLAS dgemm) with 4.7, 3.1 and 2.6 % of total time respectively > > NOTE however that gprof profiler paints to an extent different picture and > does not even mention external BLASS functions such as dtrsv and ztrsv > solvers reported as heavy users by the VerySleepy "external" profiler. > > All in all, it appears form the both profiler reports and my initial > inspection that , for start (and as I initially already intended and > suggested to), we should refactor the current heavy use of the un-productive > General Matrix copy constructor and its current reliance on > element-by-element get() function before we get into any further performance > improvements of the productive functions and external libraries. > > Best regards > > George > > ----- Original Message ----- > From: "Michel Juillard" <michel.juillard(a)ens.fr> > To: "List for Dynare developers" <dev(a)dynare.org> > Sent: Saturday, June 06, 2009 3:08 PM > Subject: Re: [DynareDev] Kalman Filter > > > >> There are tools to do profiling in C++. All we need is an standalone >> executable calling the filter. Don't loose time adding timing function >> inside the code. It may be difficult to do profiling in Windows. In that >> case, just prepare the code and we will do the profiling in Linux. >> >> Best >> >> Michel >> >> >> G. Perendia wrote: >> >>> Dear Michel >>> >>> 1) Yes, as agreed initially >>> >>> 2) these are the Matlab Dynare KF measures, mainly to show the >>> > proportion > >>> of inversion vs. pure update in Matlab KF. >>> I have not yet done fine profiling for C++, so, not much to upload >>> > either. > >>> 3) I agree.. >>> >>> Best regards >>> >>> George >>> >>> >>> ----- Original Message ----- >>> From: "Michel Juillard" <michel.juillard(a)ens.fr> >>> To: "List for Dynare developers" <dev(a)dynare.org> >>> Sent: Saturday, June 06, 2009 2:10 PM >>> Subject: Re: [DynareDev] Kalman Filter >>> >>> >>> >>> >>>> Thanks George >>>> >>>> One of the first thing that we need to establish is whether identical >>>> basic matrix operations take much longer in the C++ implementation than >>>> in Matlab and, if that it is the case, why. >>>> >>>> >>>> >>>>> 2) Indeed, and as a significant part of the overall parcel of updating >>>>> >>>>> >>> P, >>> >>> >>>>> one needs to invert the updated F too : >>>>> >>>>> 100000 loops of small model KF 4x4 F matrix inversion: iF = >>>>> >>>>> >>> inv(F); >>> >>> >>>>> Fmx_inv_time = 2.2530 >>>>> >>>>> 100000 loops of the corresponding core KF 8x8 P matrix update: >>>>> P1 = T*(P-K*P(mf,:))*transpose(T)+QQ; >>>>> Pupdt_time = 3.4450 >>>>> >>>>> (and also, 100000 loops of the preceding K = P(:,mf)*iF; >>>>> Kupdt_time = 0.5910) >>>>> >>>>> >>>>> >>>>> >>>> How do these operations compare with Matlab on your machine? >>>> >>>> >>>> >>>>> The convergence of P exploited in Matlab Dynare KFs (which does not >>>>> >>>>> >>> require >>> >>> >>>>> further update of P and K and inversion of F ), can improve greatly >>>>> performance >>>>> of KF. >>>>> >>>>> e.g.: running Matlab Dynare KF with 57x57 system matrix in 1000 loop >>>>> 1000 of usual matlabKF_time = 337.1650 >>>>> >>>>> and then using P recursively in the loop with a modified >>>>> > kalman_filter.m > >>>>> which returns P too (therefore, utilising P convergence and avoiding >>>>> > its > >>>>> update for the most of the remaining 999 loops): >>>>> 1000 of recursive: Matlab_rec_KF_time = 11.7060 >>>>> >>>>> >>>>> >>>>> 3) And, although the convergence of P in Matlab KF did not take place >>>>> >>>>> >>> for >>> >>> >>>>> the large sw_euro model which had total run much closer to the C++KF >>>>> >>>>> >>> (see 3 >>> >>> >>>>> below), as today's check show, the convergence did however take place >>>>> >>>>> >>> very >>> >>> >>>>> early in the Matlab KF running the small model I initially tested (at >>>>> >>>>> >>> step >>> >>> >>>>> t=3!!!) so it >>>>> certainly did affect and, judging from the above results, rather >>>>> >>>>> >>> greatly >>> >>> >>>>> contribute to the very much faster KF loops we experienced running >>>>> Matlab KF versus C++ KF in the initial tests with the same small >>>>> > model > >>>>> (C++ KF does yet not take >>>>> advantage of convergence and the comparative results were even)!!! >>>>> >>>>> >>>>> >>>>> >>>> OK, we forget the first comparison on the small model, because C++ and >>>> Matlab didn't use the same algorithm (no convergence monitorinig in >>>> > C++). > >>>> Matlab is still faster by 45% on the medium size model. We should focus >>>> on explaining this difference and we don't need to bring monitoring the >>>> convegence of the filter for this particular example. >>>> >>>> Could you please upload on SVN the code that you use for profiling? >>>> >>>> Best >>>> >>>> Michel >>>> >>>> >>>> >>>>> Best regards >>>>> >>>>> George >>>>> >>>>> >>>>> >>>>> >>>> _______________________________________________ >>>> Dev mailing list >>>> Dev(a)dynare.org >>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>> >>>> >>> >>> _______________________________________________ >>> Dev mailing list >>> Dev(a)dynare.org >>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>> >>> >>> >> _______________________________________________ >> Dev mailing list >> Dev(a)dynare.org >> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>

2 15

Re: [DynareDev] returning arrays from fortran function to C
by G. Perendia 18 Jun '09

18 Jun '09

PS: As an afterthought I think a more productive and optimal alternative is to write 2-3 high level f90 shell SUBROUTINES calling the QT functions and performing the several operations needed to construct the resulting Ta and TPT' and those will then also act as interface with C++. Best regards George artilogica(a)btconnect.com ----- Original Message ----- From: "G. Perendia" <george(a)perendia.orangehome.co.uk> To: "List for Dynare developers" <dev(a)dynare.org> Sent: Thursday, June 18, 2009 8:06 PM Subject: returning arrays from fortran function to C > Dear Andrea and Michel (et al) > > I have encountered an unexpected problem integrating KalmanFilter with the > f90 QT library - passing the QT result arrays back to C++. > > QT Fortran routines have been written in standard Fortran FUNCTION format, > (i.e., not SUBROUTINE), so that they are returning double or single > dimensional array (they are named by), by value ( not reference) and, as it > appears, only simple, single variables seems can be passed from Fortran > FUNCTIONs back to C++ (e.g. INT or REAL). > > On the other hand, NAG, BLAS and LAPACK routines have all been written as > Fortran SUBROUTINEs and the subroutines are the equivalent of "C" functions > returning "(void)". > > SUBROUTINE can be integrated with C more easily as they receive parameters > and return their results through the variables passed as calling parameters > by references. > > For example, dgemv.f from BLAS library gets Y by reference and returns > modified Y passed as calling parameter reference. > > SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) > .... > * Y - DOUBLE PRECISION array of DIMENSION at least > ... > * Before entry .... the incremented array Y > * must contain the vector y. On exit, Y is overwritten by the > * updated vector y. > .... > > I could not find any references on how to get arrays from Fortran FUCTION as > return value back to C - does anyone know how to do it or, if at all > possible? > > One way I can think of is less explored option of returning Fortran pointer > to the resulting array from the QT functions instead of the array by value > and I think I can work one of that out. > > However, even if there is another way to pass function result arrays back to > C++, I expect it is bound to be less efficient than passing it by reference > (or Fortran pointer), especially for larger matrices. And, if there is no > efficient alternative, I can either easily rewrite QT library as SUBROUTINE > routines instead FUCTION or try to use Fortran pointers and, if ok, I can > then also rewrite QT library to return pointers.. > > Best regards > > George >

1 0

Re: [DynareDev] Modfile en html
by Sébastien Villemot 17 Jun '09

17 Jun '09

Pablo Winant <pablo.winant(a)gmail.com> a écrit : > Je suis tout à fait d'accord pour mathml. C'est globalement mieux. Il > faut juste faire un choix ente le mode sémantique et le mode > présentation. Pour Dynare, je crois qu'il est plus facile d'écrire le > mode sémantique, qui est de plus celui indiqué pour sérialiser le modèle > (en regardant vite fait j'ai l'impression que sérialiser directement le > modfile demande des changements un peu plus intrusifs). Je peux faire un > test rapide de rendu d'un modèle en Mathml, sans toucher au code source > et voir quel rendu on a. > > jsmath est sous license apache 2 c'est écrit sur la page de Download, > mais je n'ai jamais pensé qu'on pouvait le distribuer. Je me suis juste > dit que c'était plus facile à installer et à utiliser qu'une > distribution Latex nécessaire pour compiler le fichier .tex. Et puis ça > évite dans un premier temps de réécrire la partie correspondante dans > ExprNode.cc. > > Du reste, j'ai surtout fait ça pour me faire la main, si Mathml > fonctionne bien je trouve aussi que c'est la meilleure solution. Désolé > pour les changements d'indentation, il faut que je retienne mon éditeur > de faire ça tout seul ! Je pense que ce serait intéressant d'intégrer la fonctionnalité HTML avec le support MathML. Donc si tu as le temps, n'hésite pas à l'implémenter et j'appliquerai le patch. A+ -- Sébastien Villemot

1 0

Dynare 4.0.4
by Sébastien Villemot 17 Jun '09

17 Jun '09

Hi all, We are planning to release Dynare 4.0.4 on friday, in order to have a new version for the summer school. Please merge any bugfix that you made in the trunk and which is applicable to the 4.0 branch. Best, -- Sébastien Villemot

1 0

Re: [DynareDev] Modfile en html
by Sébastien Villemot 17 Jun '09

17 Jun '09

Salut Pablo, Pablo Winant <pablo.winant(a)gmail.com> a écrit : > La dernière fois qu'on s'est vu, je t'avais parlé d'une librairie > Javascript jsmath (http://www.math.union.edu/~dpvc/jsMath/) qui permet > d'afficher des formules latex dans des pages html moyennant > l'installation des fontes latex dans le navigateur. > > En regardant comment tu as implémenté l'option pour sortir des fichiers > latex, j'ai écrit un petit patch (joint) qui ajoute la même chose pour > du html avec les commandes 'write_html_static_model' et > 'write_html_dymamic_model'. Il écrit le fichier dans le répertoire > courant et par convention, il s'attend à trouver dans le même > répertoire, un dossier jsMath (téléchargeable), et un fichier de style > 'style.css' (joint et très améliorable). Je pense que ce chemin devrait > être spécifiable par une option, mais je ne sais pas trop comment faire > pour l'instant. Merci pour le patch, l'idée est intéressante. Je me demandais pourquoi ne pas utiliser MathML? C'est a priori la façon standard d'insérer des maths dans une page web (recommandée par le W3C). Firefox le supporte directement. Pour IE il faut mettre un plugin. Dans les deux cas j'ai l'impression que l'installation est plus simple que pour jsMath. Voir: http://en.wikipedia.org/wiki/Mathml Le gars qui fait jsMath dit que MathML c'est bien, mais qu'il a fait jsMath parce que c'est plus simple à écrire pour un humain. Mais justement on génère le fichier via un programme. Donc quel avantage à utiliser jsMath qui est non-standard? Par ailleurs je n'ai pas réussi à déterminer la licence sous laquelle est distribué jsMath. Sinon un détail: le patch que tu m'as envoyé est constitué pour une bonne part de changements dans l'identation du code, changements qui cassent les conventions d'indentation C++ pour Dynare. Voir: http://www.dynare.org/DynareWiki/CodingStandards A+ -- Sébastien Villemot

1 0

Re: [DynareDev] Kalman Filter
by Michel Juillard 06 Jun '09

06 Jun '09

One more thing: what is the size of your testing model. How many variables total ? How many state variables? How many observed variables? Can you redo the experiments 3 and 4 for a model with more than 100 state variables? Probably 100 calls will be enough Best Michel G. Perendia wrote: > Dear Michel > > Results of run lapse time using etime(clock,t) are rather surprising and > very disappointing (I had to repeat some in a sheer disbelief) - far the > best is Matlab Dynare kalman_filter.m (even after applying some refactoring > of the kalman.cpp's deep inner loop, not applied to all DLL tests yet)!!!. > > Timing tests on Pentium 2.4 GHz > > 1)calling Kalman_filters.m in the 10000 loop, itself calling the optimised > version of DLL with preparation of H and Z matrices in each loop: > total_dll_time = 202.7320 > > 2)calling core optimised version of Kalman_filter_dll in the 10000 loop > after (i.e., loop without) preparation of H&Z matrices in kalman_filters.m: > core_dll_time= = 1st run: 172.5890 2nd run: 195.0500 at ~93% of CPU time > > 3)calling the Kalman_filter_loop.dll with its own inner 10000 loop from > kalman_filters.m (including the one-off call time) > innrloop_dll_time = 1st run 197.8640 2nd run: 192.5970 > 3rd run (after some refactoring of the kalman.cpp filter deep inner loop > not yet used in other DLL tests): 184.9760 > > 4)calling Dynare SVN kalman_filter.m in 10000 loop > matlabKF_time = 1st run: 48.9600 2nd run: 48.4600 > > 5) for a comparison and sanity check: calling Kalman_filters running the > debug version of DLL in the 10000 loop > dbug total_dll_time=etime(clock, t) = 642.4140 running at ~93% of CPU time > > It is possible that using proper, optimised lapack and BLAS may give a > better performance for C++. (NOTE: the DLLs used Matlab lapack and blas > libraries.) > > Shall I nevertheless continue with the C++ version of DsgeLikelihood as > planned or try to revise the performance tests first? > > Best regards > > George > artilogica(a)btconnect.com > > ----- Original Message ----- > From: "Michel Juillard" <michel.juillard(a)ens.fr> > To: "G. Perendia" <george(a)perendia.orangehome.co.uk> > Sent: Monday, June 01, 2009 10:32 AM > Subject: Re: [DynareDev] Kalman Filter > > > >> Thanks George, now it is clear. I'm also interested by a timing >> comparison. Could you time >> a) 10'000 calls to the Matlab filter routine >> b) 10'000 calls to the C++ filter routine >> c) add the 10'000 loop inside the DLL and time one call to this modified >> function. >> >> I guess after that we should pass to the coding in C++ of DsgeLikelihood.m >> >> All the best, >> >> Michel >> >> >> G. Perendia wrote: >> >>> Dear Michel >>> >>> No, I do not think there is additional problem (other than with Diffuse >>> Likelihood*.m files when if used) or that correction was not sufficient. >>> >>> Over the weekend I was just trying full 1st stage estimation with a >>> > model we > >>> use for Part Info, with and without the corrections, which just shows >>> > the > >>> scale of the effect of the log likelihood calculation error (now >>> > properly > >>> corrected for the SVN likelihood files) . As you can see in the .rtf >>> > file, > >>> in the few of the tests with SVN version (using >>> > /likelihood/kalman.filter.m) > >>> I used your new code too, but for comparing results in older versions of >>> Dynare- v4.0.2 and 4.0.3 - (using Diffuse Likelihood1.m), I used my >>> > quick > >>> fix applied to the 4.02 DiffuseLikelihod*.m that I also enclosed. >>> >>> Also, some additional tests I did over the weekend show that, after the >>> correction, Dynare now reports estimates more comparable to the new Part >>> Info method than before the correction, also proving that the correction >>> works ! >>> >>> As you suggested, I run a test of your new correction that can be >>> > compared > >>> to C++ and, for the test model I used, F converges at t=3 (start is 41), >>> > and > >>> your new kalman_filter.m code gives total >>> LIK=1640935.5855619563 >>> which is virtually the same as C++ and the quick fix for Diffuse >>> > Likelihood > >>> I tested last week (requoted from below): >>> >>> >>> >>>>>>> the reported likelihood for presample=40 from (quick >>>>>>> >>>>>>> >>> fix) Matlab KF >>> >>> >>>>>>> 1640935.5855267849 >>>>>>> which is nearly the same as that from C++ KF below: >>>>>>> 1640935.5854489324 >>>>>>> >>>>>>> >>> Let me know if you want me to run more tests or work on correcting any >>> > other > >>> Dynare Matlab code! >>> >>> Best regards >>> >>> George >>> >>> ----- Original Message ----- >>> From: "Michel Juillard" <michel.juillard(a)ens.fr> >>> To: "G. Perendia" <george(a)perendia.orangehome.co.uk> >>> Sent: Monday, June 01, 2009 9:08 AM >>> Subject: Re: [DynareDev] Kalman Filter >>> >>> >>> >>> >>>> Dear George, >>>> >>>> I'm a bit confused. We agreed that all the versions of Kalman >>>> likelihood, both in ./matlab and in ./matlab/kalman/likelihood were >>>> wrong when presample was larger than the period where the filter >>>> > reaches > >>>> its stationary state. So I don't see the point in still using this code >>>> for testing. >>>> >>>> the code in matlab/kalman/likelihood was corrected last Wednesday on >>>> > SVN > >>>> r2007 >>>> >>>> Do you have reasons to believe that the correction wasn't sufficient? >>>> > Or > >>>> is there another source of discrepancy between the Matlab code and the >>>> C++ code? >>>> Before comparing results of entire estimation procedures, we have to >>>> make sure that all the outputs from the Matlab and C++ functions with >>>> the same inputs are close enough. >>>> >>>> Note that we still have to correct DsgeLikelihood_hh.m that calls the >>>> obsolete DiffuseLikelihood*.m routines instead of the ones in >>>> kalman/likelihood >>>> >>>> Then we should just remove the DiffuseLikelihood*.m functions >>>> >>>> >>>> Best >>>> >>>> Michel >>>> >>>> >>>> G. Perendia wrote: >>>> >>>> >>>>> Dear Michel >>>>> >>>>> Problem and discrepancies with log likelihood calculation with >>>>> >>>>> >>> presampling >>> >>> >>>>> (start>1) is something I have been encountering for some time (i.e. >>>>> > few > >>>>> weeks) already , whilst doing some tests for the new partial >>>>> > information > >>>>> kalman estimation, but, expecting that the problem is in the new >>>>> > Partial > >>>>> info code, I did not get to the bottom of the cause until I >>>>> > encountered > >>>>> discrepancy between the C++ KF (modified to handle presampling >>>>> > properly) > >>> and >>> >>> >>>>> the Dynare KF too. It was that second discrepancy that led me to >>>>> > analyse > >>>>> Dynare KF handling of presampling and to find the problem there. >>>>> >>>>> I did some more preliminary tests yesterday with the older and the >>>>> >>>>> >>> latest >>> >>> >>>>> SVN versions of Dynare, without and with the correction (including >>>>> correction using the fix based around the min() function I initially >>>>> suggested and now applied as a quick fix to the older >>>>> >>>>> >>> DiffuseLikelihood1.m) >>> >>> >>>>> and the estimation results with and without correction are >>>>> > significantly > >>>>> different as can be seen from the enclosed rtf document which uses two >>>>> >>>>> >>> very >>> >>> >>>>> similar models, (i.e. one with gamma>0 the other with gamma=0), to >>>>> >>>>> >>> compare >>> >>> >>>>> the results (also very similar to the one I sent to Stephane). >>>>> >>>>> Although the initial calculations of log-likelihood ("Initial value of >>>>> >>>>> >>> the >>> >>> >>>>> log posterior (or likelihood): ") are same, (i.e. when reste is 0 or >>>>> >>>>> >>> very >>> >>> >>>>> small), it appears that the differences builds up during the >>>>> > estimation > >>> as >>> >>> >>>>> the kalman filter cov matrix start to converges earlier and earlier, >>>>> eventually before the given presample start is reached, as the later >>>>> > in > >>> the >>> >>> >>>>> ML and estimation error minimisation we are (and probably throughout >>>>> > the > >>>>> MCMC estimation process), when the parameters get close to the mode >>>>> (especially if the parameter priors were taken from previously >>>>> > estimated > >>>>> posteriors as is the case with the test models I used). >>>>> >>>>> (I did not yet run MCMC due to estimation failing with the chol() >>>>> non-positive definite error so I would need to change starting >>>>> > values.) > >>>>> Let me know if you would like me to develop and perform more >>>>> > structured > >>>>> tests on a large range of models. >>>>> >>>>> >>>>> Best regards >>>>> >>>>> George >>>>> >>>>> ----- Original Message ----- >>>>> From: "Michel Juillard" <michel.juillard(a)ens.fr> >>>>> To: "List for Dynare developers" <dev(a)dynare.org> >>>>> Cc: "G. Perendia" <george(a)perendia.orangehome.co.uk> >>>>> Sent: Wednesday, May 27, 2009 12:16 PM >>>>> Subject: Re: [DynareDev] Kalman Filter >>>>> >>>>> >>>>> >>>>> >>>>> >>>>>> In fact the discrepancy found by George occurs only for start ~= 1 >>>>>> So it is clearly related to the handling of the constants. >>>>>> >>>>>> Best >>>>>> >>>>>> Michel >>>>>> >>>>>> >>>>>> Stéphane Adjemian wrote: >>>>>> >>>>>> >>>>>> >>>>>>> George, Can you share your example. In my experience (with medium >>>>>>> scaled models) it takes much more iterations to get to the steady >>>>>>> state kalman filter... >>>>>>> >>>>>>> I didn't look at the cc codes, but one possible source of >>>>>>> > discrepancy > >>>>>>> between matlab and cc may be the steady state KF criterion. I test >>>>>>> > for > >>>>>>> the stationarity of KF considering the changes in the biggest >>>>>>> > elements > >>>>>>> of the Kalman Gain matrix, not the changes in the determinant of the >>>>>>> forecast errors covariance matrix... >>>>>>> >>>>>>> Best, >>>>>>> Stéphane. >>>>>>> >>>>>>> >>>>>>> 2009/5/27 G. Perendia <george(a)perendia.orangehome.co.uk >>>>>>> <mailto:george@perendia.orangehome.co.uk>> >>>>>>> >>>>>>> >>>>>>> 1. In my experience, the filter cov. determinants dF becomes >>>>>>> stationary after several steps (5-10), e.g. 7 >>>>>>> and number of periods the filter is stationary is >>>>>>> reste=smpl-7 >>>>>>> whilst, if start =41 (presample=40) than reste is then larger >>>>>>> > than > >>>>>>> smpl -start+1=smpl-40. >>>>>>> >>>>>>> 2. it is not problematic to add all the determinants at once in >>>>>>> the last period of the stationary filter >>>>>>> as long as we subtract from the total sum >>>>>>> >>>>>>> start-1-(smpl-reste) = 41 -1 -(smpl-smpl+7)=33 >>>>>>> determinants, or, add only >>>>>>> min(reste,(smpl-start+1)) >>>>>>> Best regards >>>>>>> >>>>>>> George >>>>>>> >>>>>>> ----- Original Message ----- >>>>>>> *From:* Stéphane Adjemian >>>>>>> > <mailto:stephane.adjemian@gmail.com> > >>>>>>> *To:* List for Dynare developers <mailto:dev@dynare.org> >>>>>>> *Cc:* G. Perendia <mailto:george@perendia.orangehome.co.uk> >>>>>>> *Sent:* Wednesday, May 27, 2009 11:04 AM >>>>>>> *Subject:* Re: [DynareDev] Kalman Filter >>>>>>> >>>>>>> Hi all, >>>>>>> >>>>>>> I agree, the matlab code is very unclear (even if I had fun >>>>>>> writting it this way ;-) and prone to errors if one uses the >>>>>>> vector lik (Marco is using it). I would rather prefer to add >>>>>>> the constants outside of the loop with a (sub)vector >>>>>>> operation, this should be more efficient. I will do it today >>>>>>> or tomorrow. >>>>>>> >>>>>>> Best, >>>>>>> Stéphane. >>>>>>> >>>>>>> 2009/5/27 Michel Juillard <michel.juillard(a)ens.fr >>>>>>> <mailto:michel.juillard@ens.fr>> >>>>>>> >>>>>>> On closer inspection, I don't think that the expression >>>>>>> pointed by George in kalman_filter.m is wrong: >>>>>>> >>>>>>> 1. reste = smpl-t or the number of periods during which >>>>>>> the filter is stationary. This shouldn't be larger than >>>>>>> T-start+1 >>>>>>> >>>>>>> 2. it is problematic (see below) but not wrong to add >>>>>>> > all > >>>>>>> the determinants at once in the last period of the >>>>>>> stationary filter >>>>>>> >>>>>>> 3. I don't think this explains the difference with the >>>>>>> > C++ > >>>>>>> version of the filter and we still have to look for it. >>>>>>> >>>>>>> 4. it remains that the current code is very unclear and >>>>>>> that if LIK is correct the vector lik doesn't have the >>>>>>> correct constants on each elements. >>>>>>> >>>>>>> 5. I would like to simplify the code and add the correct >>>>>>> constant to each element of the lik vector. It would be >>>>>>> > a > >>>>>>> little bit less efficient in Matlab than the current >>>>>>> > code, > >>>>>>> but I doubt it would be noticeable. >>>>>>> Stephane, what do you think? >>>>>>> >>>>>>> >>>>>>> Best >>>>>>> >>>>>>> Michel >>>>>>> >>>>>>> G. Perendia wrote: >>>>>>> >>>>>>> Dear Michel >>>>>>> >>>>>>> I think I found an error in Dynare Matlab >>>>>>> kalman_filter. suite of utilities >>>>>>> which affects the likelihood LIK results with >>>>>>> > start>1 > >>>>>>> (i.e. presampling>0): >>>>>>> >>>>>>> the calculation speed-up construct which relies on >>>>>>> converged covariance >>>>>>> matrix >>>>>>> >>>>>>> lik(t) = lik(t) + reste*log(dF); >>>>>>> >>>>>>> adds reste * log(dF) to the last-1 (i.e. the smpl) >>>>>>> member of lik >>>>>>> (the last, the lik(smpl+1) one contains >>>>>>> >>>>>>> >>>>>>> >>>>> smpl*pp*log(2*pi)) >>>>> >>>>> >>>>> >>>>>>> but reste is usually larger than T-start+1 so that >>>>>>> >>>>>>> LIK = >>>>>>> >>>>>>> >>>>>>> >>>>> .5*(sum(lik(start:end))-(start-1)*lik(smpl+1)/smpl) >>>>> >>>>> >>>>> >>>>>>> has much more log(dF)s added than required since >>>>>>> > they > >>>>>>> are all concentrated >>>>>>> in the last-1 (the T) member >>>>>>> >>>>>>> For example, if I change the above construct to >>>>>>> lik(t) = lik(t) + min(reste,(smpl-start+1))*log(dF); >>>>>>> >>>>>>> the reported likelihood for presample=40 from Matlab >>>>>>> >>>>>>> >>> KF >>> >>> >>>>> is >>>>> >>>>> >>>>> >>>>>>> 1640935.5855267849 >>>>>>> which is nearly the same as that from C++ KF below: >>>>>>> 1640935.5854489324 >>>>>>> >>>>>>> Shall I make changes to kalman/likelihood/ KFs and >>>>>>> upload the .m files? >>>>>>> This problem affects also the older versions of >>>>>>> DiffuseLikelihood**.m too. >>>>>>> >>>>>>> Best regards >>>>>>> >>>>>>> George >>>>>>> artilogica(a)btconnect.com >>>>>>> >>>>>>> >>>>>>> >>>>> <mailto:artilogica@btconnect.com> >>>>> >>>>> >>>>> >>>>>>> ----- Original Message ----- From: "Michel Juillard" >>>>>>> <michel.juillard(a)ens.fr >>>>>>> >>>>>>> >>> <mailto:michel.juillard@ens.fr>> >>> >>> >>>>>>> To: "G. Perendia" <george(a)perendia.orangehome.co.uk >>>>>>> <mailto:george@perendia.orangehome.co.uk>> >>>>>>> Sent: Tuesday, May 26, 2009 10:32 AM >>>>>>> Subject: Re: Kalman Filter+PS >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> Hi George, >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> Re 1) below: >>>>>>> I modified C++ KF so that it reports >>>>>>> log-likelihood for given >>>>>>> start/preampling in same/similar manner as >>>>>>> > the > >>>>>>> Matlab KFs do and I am >>>>>>> getting approximately close results, e.g. >>>>>>> >>>>>>> ll= -1640935.5854489324 for C++ and >>>>>>> (-) 1640482.4179242959 for Matlab KF (for >>>>>>> start=41, i.e. >>>>>>> >>>>>>> >>>>>>> presample=40). >>>>>>> >>>>>>> >>>>>>> whilst they appear same for presample=0 >>>>>>> (e.g.2.5906e+006), i.e. >>>>>>> -2590556.989730841 vs >>>>>>> 2590556.989778722 >>>>>>> >>>>>>> Are those results acceptably close or should >>>>>>> > I > >>>>>>> investigate further where >>>>>>> >>>>>>> >>>>>>> the >>>>>>> >>>>>>> >>>>>>> above difference may come form? >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> This indicates a problem . The difference should >>>>>>> be the same with and >>>>>>> without presample. It may come from the >>>>>>> computation of the likelihood >>>>>>> constant. This is done in a very obscure manner >>>>>>> > in > >>>>>>> Dynare Matlab. >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> _______________________________________________ >>>>>>> Dev mailing list >>>>>>> Dev(a)dynare.org <mailto:Dev@dynare.org> >>>>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> -- >>>>>>> Stéphane Adjemian >>>>>>> CEPREMAP & Université du Maine >>>>>>> >>>>>>> Tel: (33)1-43-13-62-39 >>>>>>> >>>>>>> >>>>>>> >>>>>>> > ------------------------------------------------------------------------ > >>>>>>> _______________________________________________ >>>>>>> Dev mailing list >>>>>>> Dev(a)dynare.org <mailto:Dev@dynare.org> >>>>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>>>> >>>>>>> >>>>>>> _______________________________________________ >>>>>>> Dev mailing list >>>>>>> Dev(a)dynare.org <mailto:Dev@dynare.org> >>>>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> -- >>>>>>> Stéphane Adjemian >>>>>>> CEPREMAP & Université du Maine >>>>>>> >>>>>>> Tel: (33)1-43-13-62-39 >>>>>>> >>>>>>> >>>>> ----------------------------------------------------------------------- >>>>> > - > >>>>>>> _______________________________________________ >>>>>>> Dev mailing list >>>>>>> Dev(a)dynare.org >>>>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>> _______________________________________________ >>>>>> Dev mailing list >>>>>> Dev(a)dynare.org >>>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>>> >>>>>> >>>>>> >>>>>> >>> >>> >>> >> > > > >

2 9

[Fwd: Quasi triangular matrices in Kalman filters]
by Michel Juillard 01 Jun '09

01 Jun '09

Here are the routines prepared by Andrea All the best, Michel -------- Original Message -------- Subject: Quasi triangular matrices in Kalman filters Date: Mon, 1 Jun 2009 20:47:20 +0200 From: andrea pagano <pagano.andrea(a)gmail.com> To: george(a)perendia.orangehome.co.uk CC: Michel Juillard <Michel.Juillard(a)ens.fr> Hi all I am sending you a set of Fortran routines to calculate the matricial expression in Kalman filter together with some explenations. Hope they can be a starting point in optimizing the overall computation Best Andrea -- Andrea Pagano via Veratti VARESE tel. +3903321691261 cell.+393403804397

1 0