bugGNU Octave - Bugs: bug #52569, unexptected error occured when...

 
 

bug #52569: unexptected error occured when using dynare with octave

Submitter:  None
Submitted:  Fri 01 Dec 2017 01:33:22 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Klára Major Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 4.2.1
Operating System:  * Microsoft Windows Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 03 Mar 2019 03:00:17 PM UTC, comment #8: 

Excellent.  Marking as fixed and closing report.

Rik <rik5>
Group administrator
Sun 03 Mar 2019 08:36:35 AM UTC, comment #7: 

Indeed the problem does not persist in Octave 5.1

Johannes Pfeifer <jpfeifer>
Sun 03 Mar 2019 05:20:50 AM UTC, comment #6: 

Version 4.2.1 is deprecated and no longer supported.  Please try with the latest stable release which is 5.1.0.  The problem is likely to have been fixed.

Rik <rik5>
Group administrator
Tue 20 Mar 2018 07:51:54 AM UTC, comment #5: 

Someone might want to test this again with the upcoming Octave 4.4 on Windows. The error handling from Fortran routines was especially troublesome on Windows, but Octave's internal error handling has been completely redone for 4.4. I don't have access to Windows to test this, but it may be an improved error reporting experience.

Mike Miller <mtmiller>
Group Member
Fri 01 Dec 2017 08:25:08 PM UTC, comment #4: 

Yes, the error handling is a bit clumsy.  That's due to the fact that the core numerical routine is in FORTRAN (not my emphasis, but stands for FORMula TRANslation) so the only convenient error messages are from FORTRAN's messaging.

I can make a guess as to what the issue is with this routine.  Note the following about this part of the algorithm:


      DO 20 K=1,200
        FK = K
        T = (A+FK)*X*(1.D0+R)
        R = T/((AX+FK)*(A1X+FK)-T)


It seems that if both A and X are large, R doesn't become small (i.e., such that P = R*P starts to diminish) until K gets on the order of A and X.  So, obviously if we choose A and X >> 200 it could be that not much happens to the final value until K approaches 200.  In other words, if I were to investigate I'd look at whether the issue here isn't so much "convergence" the way we think, but more a case of not choosing the range of K suitable to the values of X and A.  Stated differently, it may be that for a given X and A there is a value of K which most contributes to the final result S = S + P and that value of K (call it K_0) may be larger than 200.  So one wonders if the problem should be transformed such that loop-indexing is more like

K_0
K_0 + 1
K_0 - 1
K_0 + 2
K_0 - 2

etc.  Something like that might always converge within a reasonable number of iterations.

However, I don't want to go into that much depth of investigation at the moment.

Dan Sebald <sebald>
Fri 01 Dec 2017 07:45:36 PM UTC, comment #3: 

Dear Dan,
thanks for looking into this. I guess the expectation was that Octave is able to solve this numerical problem the way Matlab does. Because Matlab provides a solution and Octave crashes during execution, we flagged this as a bug in Dynare.

It's difficult to say what the expectation is, given that we are talking about a computational limitation. But I would probably have expected a better error handling with a more informative error message instead of a cryptic crash. One may also debate whether one should return a warning message that the problem could not be solved and return a NaN instead of throwing an exception.

As with respect to the original user, there is kind of a workaround. Using


plot_priors=0,prior_trunc=0

in the estimation options of Dynare will prevent calling the problematic code. That will come at the cost of not being able to plot the prior distribution, but allows running the estimation. I will adjust Dynare's error message accordingly.

Johannes Pfeifer <jpfeifer>
Fri 01 Dec 2017 06:43:10 PM UTC, comment #2: 

Thanks for the links.

I've looked into the code a bit, just to get a handle on what the issue is.  SLATEC D9LGIT is a very short routine, so not difficult to comprehend.  SLATEC is a collection of lesser libraries, of which d9lgit belongs to a library called The Fullerton Function Library FN, which has been updated and ported to C/C++ in some way but I didn't look at much of that:

https://people.sc.fsu.edu/~jburkardt/f77_src/fn/fn.html
https://people.sc.fsu.edu/~jburkardt/cpp_src/fn/fn.html
http://www.netlib.org/fn/d9lgit.f

The last reference is the actual, more up-to-date code.  However, it isn't much different, aside from using the double versions of some of the routines.  I'm attaching a diff file that makes Octave's d9lgit.f in line with the one in the link above.  (Don't apply this to the repository, though.  I don't think we want to start freely modifying third-party-sourced code.)  I don't see any difference in convergence behavior when changing the ABS to DABS.

The issue is that poisscdf (A,A+E) doesn't converge well when A grows large and E shrinks, e.g., poisscdf(1000,1001) versus poisscdf(1000,1002).  I don't have time to think of the mathematical ramifications of this property (e.g., does this algorithm outright lose accuracy for the mathematical function as A becomes large?).

That being the case, it appears that simply bumping up the allowable number of iterations at least avoids the FORTRAN error.  That is, I made this change:


diff --git a/liboctave/external/slatec-fn/d9lgit.f b/liboctave/external/slatec-f
--- a/liboctave/external/slatec-fn/d9lgit.f
+++ b/liboctave/external/slatec-fn/d9lgit.f
@@ -46,7 +46,7 @@ C
       R = 0.D0
       P = 1.D0
       S = P
-      DO 20 K=1,200
+      DO 20 K=1,500
         FK = K
         T = (A+FK)*X*(1.D0+R)
         R = T/((AX+FK)*(A1X+FK)-T)
@@ -55,7 +55,7 @@ C
         IF (ABS(P).LT.EPS*S) GO TO 30
  20   CONTINUE
       CALL XERMSG ('SLATEC', 'D9LGIT',
-     +   'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2)
+     +   'NO CONVERGENCE IN 500 TERMS OF CONTINUED FRACTION', 3, 2)
 C
  30   HSTAR = 1.0D0 - X*S/A1X
       IF (HSTAR .LT. SQEPS) CALL XERMSG ('SLATEC', 'D9LGIT',
(END)


and can then get results for larger A.  For example:


octave:5> poisscdf (4000,4000)
ans =  0.50421
octave:6> poisscdf (5000,5000)
 ***MESSAGE FROM ROUTINE D9LGIT IN LIBRARY SLATEC.
 ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  NO CONVERGENCE IN 500 TERMS OF CONTINUED FRACTION
 *  ERROR NUMBER = 3
 *
 ***END OF MESSAGE

 ***JOB ABORT DUE TO FATAL ERROR.
0          ERROR MESSAGE SUMMARY
 LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
 SLATEC     D9LGIT     NO CONVERGENCE IN 50         3         2         2

error: gammainc: exception encountered in Fortran subroutine xgammainc_
error: called from
    poisscdf at line 60 column 12
octave:6>


So, obviously this algorithm is going to fall apart at some point.  Perhaps the original submitter could give us some insight as to reasonable expectations with regard to applications.

This "error: Due to a bug in Octave" isn't quite a bug, it would seem.  It's more like a limitation on numerical methods for approximating the function which perhaps we could stretch out but never totally eliminate.

(file #42546)

Dan Sebald <sebald>
Fri 01 Dec 2017 02:56:43 PM UTC, comment #1: 

This is an old won't fix issue in Octave. The relevant items are the previous bug reports https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=493869 and
https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=418158

Unless the Octave team reconsiders their stance, you should consider this as a limitation of Octave.

Johannes Pfeifer <jpfeifer>
Fri 01 Dec 2017 01:33:22 PM UTC, original submission:  

Hello,

I was using Dynare 4.5.3 for bayesian estimation of dsge models, previously mostly with matlab. Now I rerun my code now in Octave 4.2.1, officially stated at dynare website that comaptible with dynare 4.5.3. Unfortunately, the same code did not work, now I got the following error messages:

                                                                                                                                         
 ***MESSAGE FROM ROUTINE D9LGIT IN LIBRARY SLATEC.
 ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED

  •  NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION
  •  ERROR NUMBER = 3

 *
 ***END OF MESSAGE

 ***JOB ABORT DUE TO FATAL ERROR.

  1.          ERROR MESSAGE SUMMARY

 LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
 SLATEC     D9LGIT     NO CONVERGENCE IN 20         3         2         9

error: Due to a bug in Octave, you must choose other values for mean and/or variance of your prior on theta, or use another shape


Finally, by modifying the variance of theta parameter in the prior distribution, the code worked. But I just do not understand what makes the difference? Should I turn to Octave developers group or to Dynare? May I expect that it will change in the upcoming releases?

Thank you very much for your attention.

Best regards,

Klára Major

ps. I attach the files if you want to check it.


Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #42544:  dsge_est.mod added by None (1KiB - audio/wav)
file #42545:  usa_fred.mat added by None (6KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by sebald (Updated the item)
  • -email is unavailable- added by jpfeifer (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-03-03 rik5 StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2019-03-03 mtmiller Carbon-CopyRemoved 80942 -
    2019-03-03 rik5 StatusNone Need Info
    2018-03-20 mtmiller CategoryNone Libraries
        Item GroupMatlab Compatibility Unexpected Error or Warning
    2017-12-01 sebald Attached File- Added octave-d9lgit_convergence-djs2017dec01.diff, #42546
    2017-12-01 None Attached File- Added dsge_est.mod, #42544
        Attached File- Added usa_fred.mat, #42545

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code