v4 : adding Marco's sensitivity routines'

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@609 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2006-01-19 06:37:21 +00:00
parent 6e23ddbf6d
commit ec1f4ab8c3
5 changed files with 877 additions and 0 deletions

512
matlab/LPTAU.m Normal file
View File

@ -0,0 +1,512 @@
function VECTOR = LPTAU(I, N)
%
% I.M. SOBOL', V.I. TURCHANINOV, Yu.L. LEVITAN, B.V. SHUKHMAN
% KELDYSH INSTITUTE OF APPLIED MATHEMATICS
% RUSSIAN ACADEMY OF SCIENCES
%
% QUASIRANDOM SEQUENCE GENERATORS
% -------------------------------
%
% 28.11.1991
%
% NOTE TO THE USER BY the NEA Data Bank:
% This quasi random number generator has been made available to
% you on condition that its identity is preserved when used
% in computer programs. If its use leads to scientific publication
% of results you should cite it in the references, in addition
% no commercial use should be made unless agreed upon with the
% main author (Prof. I.M. Sobol')
%
% ABSTRACT
% ........
%
% POINTS BELONGING TO LP-TAU SEQUENCES UNIFORMLY DISTRIBUTED IN THE
% N-DIMENSIONAL UNIT CUBE ARE OFTEN USED IN NUMERICAL MATHEMATICS:
%
% - AS NODES FOR MULTIDIMENSIONAL INTEGRATION;
% - AS SEARCHING POINTS IN GLOBAL OPTIMIZATION;
% - AS TRIAL POINTS IN MULTI-CRITERIA DECISION MAKING;
% - AS QUASIRANDOM POINTS FOR QUASI-MONTECARLO ALGORITHMS;
% - ETC.
%
% THIS SUBROUTINE CONTAINS THE ALGORITHM FOR FAST GENERATION OF
% LP-TAU SEQUENCES THAT ARE SUITABLE FOR MULTI-PROCESSOR COMPUTATIONS.
% THE DIMENSIONS N.LE.51, THE NUMBER OF POINTS N.LT.2**30.
% THE PROGRAMMING LANGUAGE IS FORTRAN-77. THIS SUBROUTINE IS AVAILABLE
% ALSO IN %-LANGUAGE.
% THE REPORT DESCRIBING THE ALGORITHM CONTAINS THE DESCRIPTION OF THE
% ALGORITHM AND CERTAIN IMPORTANT PROPERTIES OF LP-TAU SEQUENCES AND
% THEIR GENERALIZATIONS ARE DISCUSSED.
%
% REFERENCE:
% I.M. SOBOL', V.I. TURCHANINOV, Yu.L. LEVITAN, B.V. SHUKHMAN
% KELDYSH INSTITUTE OF APPLIED MATHEMATICS
% RUSSIAN ACADEMY OF SCIENCES
%
% QUASIRANDOM SEQUENCE GENERATORS
% MOSCOW 1992, IPM ZAK. NO.30 (100 COPIES)
%
% ------------------------------------------------------------------------
%
% INPUT PARAMETERS:
%
% I - NUMBER OF THE POINT (I=(0,2**30-1)),
% N - DIMENSION OF THE POINT (0<N<52);
%
% OUTPUT PARAMETER:
%
% VECTOR(N) - N-VECTOR CONTAINING THE CARTESIAN CO-ORDINATES OF
% THE I-TH POINT.
%
%
% TO CALL THE SUBROUTINE WRITE:
%
% CALL LPTAU(I,N,VECTOR)
% WHERE I, N: INTEGER CAPABLE OF STORING 2**30 (INTEGER*4 ON IBM
% OR OTHER 32 BIT/WORD MACHINES)
% VECTOR: DOUBLE PRECISION ARRAY WHOSE LENGTH < 52.
%
% INTEGER QP
% QP = QUANTITY POWER
%
% PARAMETER (MAXDIM=51, QP=30, MAXNUM=2**30-1)
MAXDIM=51; QP=30; MAXNUM=2^30-1;
%
% THE DIMENSION OF THE POINT CANNOT EXCEED MAXDIM
% THE TOTAL NUMBER OF GENERATED POINTS CANNOT EXCEED 2**QP
% MAXNUM=2**30-1 // 1073741823
%
% DOUBLE PRECISION VECTOR(N)
% INTEGER I,N
%
% INTEGER PRVNUM,PRVDIM
% INTEGER DIRECT(MAXDIM,QP), MASKV(MAXDIM)
% DOUBLE PRECISION SCALE
%
% Translated into MATLAB by M. Ratto
VECTOR=zeros(1,N);
SCALE =9.31322574615478516E-10;
persistent PRVNUM PRVDIM MASKV DIRECT
if isempty(PRVNUM), PRVNUM=-2; end,
if isempty(PRVDIM), PRVDIM=0; end,
if isempty(DIRECT), ...
DIRECT(1:MAXDIM,1)=[
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912, 536870912, 536870912, 536870912, 536870912, ...
536870912]';
DIRECT(1:MAXDIM,2)=[
805306368, 268435456, 805306368, 268435456, 805306368, ...
268435456, 805306368, 268435456, 268435456, 805306368, ...
268435456, 805306368, 268435456, 805306368, 268435456, ...
805306368, 805306368, 268435456, 805306368, 268435456, ...
805306368, 268435456, 805306368, 268435456, 268435456, ...
805306368, 268435456, 805306368, 268435456, 805306368, ...
268435456, 805306368, 805306368, 268435456, 805306368, ...
268435456, 805306368, 268435456, 805306368, 268435456, ...
268435456, 805306368, 268435456, 805306368, 268435456, ...
805306368, 268435456, 805306368, 805306368, 268435456, ...
805306368]';
DIRECT(1:MAXDIM,3)=[
939524096, 939524096, 134217728, 671088640, 402653184, ...
402653184, 671088640, 134217728, 671088640, 402653184, ...
939524096, 134217728, 671088640, 939524096, 134217728, ...
671088640, 134217728, 939524096, 939524096, 402653184, ...
402653184, 134217728, 671088640, 402653184, 671088640, ...
402653184, 402653184, 671088640, 134217728, 134217728, ...
939524096, 939524096, 939524096, 939524096, 134217728, ...
671088640, 402653184, 402653184, 671088640, 134217728, ...
671088640, 402653184, 939524096, 134217728, 671088640, ...
939524096, 134217728, 671088640, 134217728, 939524096, ...
939524096]';
DIRECT(1:MAXDIM,4)=[
1006632960, 67108864, 603979776, 1006632960, 335544320, ...
469762048, 872415232, 738197504, 469762048, 872415232, ...
1006632960, 67108864, 872415232, 469762048, 469762048, ...
469762048, 1006632960, 335544320, 872415232, 603979776, ...
201326592, 67108864, 335544320, 67108864, 201326592, ...
738197504, 1006632960, 738197504, 603979776, 335544320, ...
738197504, 67108864, 1006632960, 67108864, 603979776, ...
1006632960, 335544320, 469762048, 872415232, 738197504, ...
469762048, 872415232, 1006632960, 67108864, 872415232, ...
469762048, 469762048, 469762048, 1006632960, 335544320, ...
872415232]';
DIRECT(1:MAXDIM,5)=[
1040187392, 637534208, 1040187392, 838860800, 167772160, ...
234881024, 167772160, 1040187392, 436207616, 33554432, ...
570425344, 1040187392, 503316480, 838860800, 973078528, ...
167772160, 234881024, 436207616, 771751936, 100663296, ...
234881024, 905969664, 771751936, 33554432, 838860800, ...
973078528, 905969664, 33554432, 301989888, 838860800, ...
100663296, 234881024, 1040187392, 637534208, 1040187392, ...
838860800, 167772160, 234881024, 167772160, 1040187392, ...
436207616, 33554432, 570425344, 1040187392, 503316480, ...
838860800, 973078528, 167772160, 234881024, 436207616, ...
771751936]';
DIRECT(1:MAXDIM,6)=[
1056964608, 352321536, 50331648, 419430400, 956301312, ...
889192448, 620756992, 117440512, 956301312, 922746880, ...
822083584, 218103808, 587202560, 385875968, 452984832, ...
218103808, 184549376, 285212672, 150994944, 956301312, ...
352321536, 654311424, 553648128, 352321536, 788529152, ...
956301312, 587202560, 251658240, 218103808, 721420288, ...
251658240, 1056964608, 520093696, 889192448, 587202560, ...
956301312, 419430400, 352321536, 83886080, 654311424, ...
419430400, 385875968, 285212672, 754974720, 50331648, ...
922746880, 989855744, 754974720, 721420288, 822083584, ...
687865856]';
DIRECT(1:MAXDIM,7)=[
1065353216, 1065353216, 578813952, 25165824, 125829120, ...
964689920, 327155712, 310378496, 360710144, 360710144, ...
159383552, 444596224, 947912704, 662700032, 444596224, ...
528482304, 578813952, 578813952, 75497472, 1065353216, ...
92274688, 511705088, 897581056, 847249408, 662700032, ...
931135488, 411041792, 713031680, 696254464, 528482304, ...
209715200, 578813952, 1065353216, 1065353216, 578813952, ...
25165824, 125829120, 964689920, 327155712, 310378496, ...
360710144, 360710144, 159383552, 444596224, 947912704, ...
662700032, 444596224, 528482304, 578813952, 578813952, ...
75497472]';
DIRECT(1:MAXDIM,8)=[
1069547520, 4194304, 826277888, 624951296, 616562688, ...
801112064, 952107008, 406847488, 398458880, 331350016, ...
356515840, 46137344, 1010827264, 1069547520, 734003200, ...
113246208, 490733568, 633339904, 910163968, 801112064, ...
893386752, 851443712, 801112064, 918552576, 742391808, ...
499122176, 62914560, 264241152, 708837376, 717225984, ...
759169024, 499122176, 272629760, 423624704, 155189248, ...
239075328, 205520896, 943718400, 373293056, 457179136, ...
406847488, 71303168, 473956352, 54525952, 624951296, ...
104857600, 1019215872, 901775360, 213909504, 281018368, ...
851443712]';
DIRECT(1:MAXDIM,9)=[
1071644672, 543162368, 190840832, 329252864, 853540864, ...
132120576, 778043392, 73400320, 178257920, 522190848, ...
639631360, 534773760, 991952896, 333447168, 48234496, ...
1059061760, 761266176, 673185792, 220200960, 396361728, ...
362807296, 815792128, 819986432, 346030080, 39845888, ...
752877568, 387973120, 643825664, 291504128, 274726912, ...
568328192, 526385152, 673185792, 98566144, 396361728, ...
727711744, 1042284544, 106954752, 299892736, 912261120, ...
44040192, 895483904, 333447168, 551550976, 467664896, ...
618659840, 606076928, 274726912, 245366784, 446693376, ...
421527552]';
DIRECT(1:MAXDIM,10)=[
1072693248, 273678336, 644874240, 753926144, 495976448, ...
869269504, 355467264, 57671680, 816840704, 961544192, ...
804257792, 495976448, 347078656, 426770432, 1066401792, ...
372244480, 84934656, 208666624, 313524224, 598736896, ...
487587840, 965738496, 1011875840, 296747008, 393216000, ...
523239424, 720371712, 823132160, 128974848, 407896064, ...
747634688, 850395136, 873463808, 504365056, 481296384, ...
686817280, 592445440, 995098624, 498073600, 969932800, ...
586153984, 1039138816, 814743552, 523239424, 294649856, ...
305135616, 506462208, 11534336, 449839104, 619708416, ...
479199232]';
DIRECT(1:MAXDIM,11)=[
1073217536, 947388416, 1070071808, 977797120, 365428736, ...
702021632, 461897728, 829947904, 425197568, 634912768, ...
437780480, 582483968, 792199168, 315097088, 611844096, ...
667418624, 166199296, 513277952, 187170816, 1036517376, ...
25690112, 201850880, 443023360, 990380032, 63438848, ...
211288064, 983040000, 1069023232, 421003264, 742916096, ...
487063552, 363331584, 973602816, 286785536, 171442176, ...
669515776, 110624768, 383254528, 289931264, 352845824, ...
878182400, 655884288, 836239360, 765984768, 549978112, ...
655884288, 85458944, 591921152, 563609600, 277348352, ...
919076864]';
DIRECT(1:MAXDIM,12)=[
1073479680, 71565312, 2359296, 891551744, 158597120, ...
383516672, 1019478016, 947126272, 621019136, 714866688, ...
738459648, 265027584, 468975616, 131858432, 504627200, ...
581173248, 266600448, 865861632, 658243584, 546045952, ...
521404416, 304873472, 1060896768, 163840000, 305922048, ...
257163264, 50069504, 773062656, 59506688, 779354112, ...
165937152, 587988992, 486801408, 160694272, 90439680, ...
423362560, 536608768, 614203392, 56885248, 999030784, ...
10747904, 764674048, 25952256, 989069312, 352583680, ...
799801344, 261357568, 873201664, 40108032, 769392640, ...
254541824]';
DIRECT(1:MAXDIM,13)=[
1073610752, 644218880, 538836992, 455475200, 1062600704, ...
139329536, 205651968, 905052160, 797048832, 452329472, ...
973471744, 627703808, 614072320, 803078144, 637403136, ...
835059712, 949878784, 662044672, 767950848, 426901504, ...
448659456, 23986176, 1016201216, 524943360, 525991936, ...
618790912, 781058048, 761659392, 458096640, 226361344, ...
950665216, 952500224, 516030464, 337510400, 496107520, ...
830865408, 944111616, 636354560, 978452480, 921567232, ...
533594112, 7471104, 678035456, 471203840, 1065746432, ...
575275008, 996540416, 909246464, 879362048, 637927424, ...
25821184]';
DIRECT(1:MAXDIM,14)=[
1073676288, 357892096, 808648704, 2424832, 2555904, ...
624230400, 69271552, 456851456, 1052966912, 600637440, ...
487260160, 794624000, 386727936, 467599360, 798031872, ...
630652928, 340983808, 493944832, 37945344, 264175616, ...
263520256, 833421312, 235077632, 464846848, 534839296, ...
992411648, 10813440, 367067136, 116457472, 115015680, ...
928710656, 619773952, 813760512, 1043398656, 967770112, ...
912850944, 72155136, 1009057792, 668532736, 462356480, ...
267321344, 795803648, 635764736, 574160896, 1003421696, ...
181075968, 56688640, 388562944, 190906368, 657915904, ...
474939392]';
DIRECT(1:MAXDIM,15)=[
1073709056, 1073709056, 137003008, 547782656, 545095680, ...
26836992, 34701312, 354385920, 925663232, 656965632, ...
327581696, 894795776, 110067712, 1038057472, 209354752, ...
596541440, 42631168, 471433216, 52527104, 666861568, ...
706707456, 674070528, 824410112, 305496064, 136282112, ...
847740928, 531464192, 222920704, 379289600, 507740160, ...
11894784, 1053392896, 129990656, 557547520, 666468352, ...
1061912576, 576684032, 1041334272, 380469248, 114196480, ...
133070848, 517046272, 129990656, 790396928, 563773440, ...
388333568, 661749760, 446791680, 737378304, 229998592, ...
348225536]';
DIRECT(1:MAXDIM,16)=[
1073725440, 16384, 605372416, 275234816, 817971200, ...
603963392, 555335680, 721534976, 997801984, 1028767744, ...
407060480, 375275520, 256688128, 1021165568, 303349760, ...
1022476288, 234143744, 106708992, 732971008, 733954048, ...
789889024, 879575040, 764657664, 762658816, 1010843648, ...
941080576, 827932672, 98942976, 1051738112, 624934912, ...
993280000, 134070272, 201375744, 567558144, 882163712, ...
649084928, 356564992, 489439232, 637091840, 60637184, ...
199278592, 815677440, 927678464, 94519296, 419184640, ...
933838848, 426655744, 911130624, 171393024, 561332224, ...
471613440]';
DIRECT(1:MAXDIM,17)=[
1073733632, 536895488, 1043685376, 679944192, 417505280, ...
301981696, 832561152, 210542592, 167501824, 1071341568, ...
229302272, 970661888, 732176384, 576659456, 402464768, ...
451584000, 368467968, 928260096, 933847040, 29319168, ...
582934528, 772612096, 330014720, 647323648, 174071808, ...
1008689152, 295919616, 353869824, 177774592, 580198400, ...
381837312, 638574592, 637558784, 679370752, 504012800, ...
747118592, 429973504, 1032609792, 932667392, 583360512, ...
969498624, 1056333824, 660955136, 247488512, 153509888, ...
242180096, 205840384, 797499392, 824565760, 234348544, ...
842326016]';
DIRECT(1:MAXDIM,18)=[
1073737728, 268455936, 52785152, 1020628992, 345018368, ...
452972544, 704442368, 255987712, 750759936, 697692160, ...
196677632, 764604416, 485625856, 522022912, 680620032, ...
362270720, 838103040, 83972096, 629133312, 46108672, ...
867561472, 725422080, 184504320, 751112192, 191918080, ...
306425856, 507310080, 30453760, 281858048, 604000256, ...
208662528, 319557632, 318779392, 476139520, 863719424, ...
567062528, 521179136, 712790016, 610299904, 293687296, ...
1023086592, 549089280, 1065242624, 707751936, 363024384, ...
16674816, 197136384, 1037561856, 195112960, 372707328, ...
992751616]';
DIRECT(1:MAXDIM,19)=[
1073739776, 939554816, 580732928, 854333440, 172619776, ...
511694848, 936142848, 518199296, 593348608, 225527808, ...
900982784, 180279296, 168904704, 62814208, 754485248, ...
730691584, 1005996032, 411174912, 249866240, 641669120, ...
1008719872, 749066240, 860993536, 94177280, 432564224, ...
226355200, 925784064, 995657728, 967731200, 436226048, ...
913799168, 549894144, 964696064, 843315200, 445863936, ...
1047422976, 548947968, 492066816, 953870336, 1002653696, ...
861440000, 385636352, 325253120, 187353088, 653584384, ...
1008269312, 748693504, 1013016576, 55814144, 255170560, ...
260708352]';
DIRECT(1:MAXDIM,20)=[
1073740800, 67126272, 829514752, 423777280, 968297472, ...
205511680, 147076096, 926669824, 202300416, 118395904, ...
381332480, 1002738688, 743042048, 292551680, 584567808, ...
284339200, 183936000, 616762368, 435221504, 159376384, ...
907322368, 595696640, 247497728, 553735168, 826051584, ...
564454400, 446024704, 214236160, 33661952, 251685888, ...
660327424, 284244992, 859868160, 722502656, 622844928, ...
324342784, 682374144, 400579584, 405353472, 605187072, ...
840682496, 212956160, 157891584, 193201152, 437990400, ...
573578240, 368053248, 580197376, 937905152, 565527552, ...
89064448]';
DIRECT(1:MAXDIM,21)=[
1073741312, 637560320, 189496832, 27474432, 129338880, ...
908054016, 641870336, 186375680, 302677504, 763663872, ...
103878144, 325187072, 858254848, 922041856, 261924352, ...
954978816, 292822528, 849512960, 210311680, 933232128, ...
691981824, 155417088, 627070464, 416795136, 182081024, ...
513433088, 848658944, 515770880, 627273216, 629169664, ...
414566912, 147450368, 698353152, 244844032, 226578944, ...
1020087808, 886978048, 389697024, 1007004160, 839646720, ...
621924864, 549962240, 609583616, 735976960, 87342592, ...
1058542080, 163066368, 307997184, 876471808, 794280448, ...
675386880]';
DIRECT(1:MAXDIM,22)=[
1073741568, 352343296, 644236032, 636735232, 615860480, ...
959444224, 287380736, 1007410432, 890187008, 399480576, ...
520092928, 643311360, 816901376, 695310080, 1019229440, ...
77034240, 733295872, 1035127552, 986582784, 332381952, ...
334852352, 364956416, 596672256, 800381696, 480316672, ...
574863104, 647347968, 702910208, 499965184, 364968704, ...
120862976, 1023256320, 995114240, 13951232, 32520448, ...
702127360, 45176064, 444945664, 237860096, 152839936, ...
530633984, 429135616, 267272448, 884808960, 933712640, ...
61605632, 174335744, 564911360, 302327552, 650589440, ...
450649344]';
DIRECT(1:MAXDIM,23)=[
1073741696, 1065385856, 1073734272, 331949184, 842310784, ...
799537536, 965852032, 369351808, 662886016, 86119808, ...
865109888, 299633792, 422735488, 181087360, 174252416, ...
1041212544, 840196224, 750314368, 391053440, 903306880, ...
742365312, 236995200, 42492800, 946000512, 771692416, ...
897405824, 613803136, 924258688, 808338304, 1038125440, ...
683814272, 177186176, 766008960, 704549248, 194555008, ...
306383744, 496592512, 416020864, 655186816, 1032204928, ...
694773632, 577910144, 45797760, 910332544, 536014976, ...
675946368, 987635840, 788223872, 353993856, 96313472, ...
85248640]';
DIRECT(1:MAXDIM,24)=[
1073741760, 4210752, 12608, 744788544, 494377792, ...
115601344, 769248448, 990895808, 851706304, 979326784, ...
692061120, 429015104, 217132864, 736067008, 55694400, ...
456152640, 631601984, 787264192, 898599104, 478383808, ...
507774272, 458270272, 392995136, 872482496, 124824768, ...
1034161344, 362141632, 1053833280, 943810496, 428920128, ...
795835200, 835462848, 961843520, 606198080, 785652672, ...
154954432, 491617344, 297351232, 580735552, 634587968, ...
185277760, 141505600, 417673280, 106907456, 395575616, ...
566958656, 352651200, 415242688, 26635200, 1030317504, ...
247212992]';
DIRECT(1:MAXDIM,25)=[
1073741792, 543187040, 536882016, 981766752, 357858144, ...
810665952, 369083488, 613015520, 86115232, 845149216, ...
606076960, 1018241504, 35245536, 635403744, 236633696, ...
407451232, 427671904, 460141792, 116344544, 990246688, ...
1024892576, 883429408, 150655392, 173476896, 197666464, ...
1016740448, 605376608, 970487008, 1006881248, 598265632, ...
1022861728, 489055392, 216680608, 371439072, 53140576, ...
965114400, 98534112, 209692256, 264598496, 321437280, ...
545303840, 324119904, 876587488, 778239712, 802033120, ...
44406496, 665829024, 803213920, 309742112, 735181344, ...
1036125600]';
DIRECT(1:MAXDIM,26)=[
1073741808, 273698896, 805312112, 903123536, 153504624, ...
689632208, 432848944, 859888752, 510853776, 240805744, ...
250610192, 852489168, 139460144, 283082224, 222702928, ...
342181488, 922872432, 187528656, 360637424, 814514736, ...
301037840, 800872624, 272595184, 158627600, 238839088, ...
927181136, 710248688, 788854608, 1048376560, 484709072, ...
85709008, 1065469744, 429237328, 490778384, 848024144, ...
443114288, 687907504, 474573648, 45706416, 681465296, ...
805303216, 14567920, 369839504, 440402480, 797652624, ...
115959088, 929784560, 91226544, 205943056, 288358704, ...
950217296]';
DIRECT(1:MAXDIM,27)=[
1073741816, 947419256, 134232008, 460100200, 1073685336, ...
398354904, 1069834248, 686054696, 108299224, 273898840, ...
731382552, 938170664, 86812552, 677346808, 602208712, ...
652635752, 209797064, 323968376, 384078968, 561178056, ...
923399256, 996597176, 942777416, 885167400, 89791048, ...
956122504, 87393464, 987661336, 993412952, 827749496, ...
903211272, 33649816, 594875176, 65052056, 822835240, ...
625704888, 1065374584, 612232920, 536299720, 1046858504, ...
939518840, 160843032, 654656088, 744496936, 872197704, ...
219111576, 829112632, 455614152, 1064192952, 313010856, ...
820754920]';
DIRECT(1:MAXDIM,28)=[
1073741820, 71582788, 603989028, 37500, 120660, ...
182195932, 213921796, 473570692, 41763004, 156352828, ...
88343188, 698012780, 666108868, 337608156, 1054329884, ...
799472220, 641885244, 392104980, 502284340, 1002405628, ...
249907140, 75586228, 206587660, 565275348, 1021426548, ...
425987996, 598058580, 401940420, 919427316, 822049324, ...
115655868, 759299588, 562394644, 990942164, 1003212268, ...
598750292, 675334852, 675293340, 445665316, 903023756, ...
872412692, 172640916, 1051615436, 91229620, 94586692, ...
344873452, 7029156, 683421900, 434258804, 955002092, ...
436424380]';
DIRECT(1:MAXDIM,29)=[
1073741822, 644245094, 1040195102, 537039474, 537089354, ...
642656334, 73402402, 664239174, 1014620426, 500917234, ...
1042677826, 347788318, 1069154738, 373259762, 362587674, ...
239395914, 132283950, 903121466, 445210638, 968851198, ...
230153254, 935549622, 308815630, 861891286, 496995094, ...
670740382, 657311830, 573565382, 248331478, 1064650798, ...
338312798, 669059078, 1065190902, 379125622, 111897166, ...
520650422, 740300390, 1046483950, 66254898, 992513134, ...
234871606, 610553330, 450553866, 566758134, 787800806, ...
1071709866, 971732606, 528102354, 790919122, 917384366, ...
656244162]';
DIRECT(1:MAXDIM,30)=[
1073741823, 357913941, 50344755, 268538457, 805540473, ...
3487029, 3146769, 592038967, 729591963, 781578389, ...
66781679, 113956361, 331153483, 967802327, 935343371, ...
324351309, 609305915, 818137857, 131059769, 918519549, ...
827341719, 922770101, 456972047, 363883221, 1057014661, ...
900956063, 580478293, 520383687, 315470533, 227601711, ...
169598765, 909227271, 796983815, 349496709, 393974911, ...
378320037, 777004343, 927999269, 616377385, 345930959, ...
989849787, 627400495, 892675125, 260314121, 1041964063, ...
531367163, 195776757, 237309785, 949187605, 719002587, ...
495745187]';
end
%
% TRAP FOR A WRONG SUBROUTINE CALL
%
if ((I<0) | (N<1) | (I>MAXNUM) | (N>MAXDIM)),
disp('LP-TAU CALL FAILED')
disp(' PRESS <ENTER> TO EXIT LPTAU')
pause
return
end
if ((PRVNUM+1==I) & (N<=PRVDIM)),
%
% RECURRENT GENERATION OF THE POINT
%
%
% SEARCH POSITION OF THE RIGHTMOST "1"
% IN THE BINARY REPRESENTATION OF I
%
L=0;
POS=0;
while L<QP & POS==0,
L=L+1;
POS=bitand(bitshift(I,-(L-1)),1);
end
%
% RIGHTMOST POSITION IS L
%
for J=1:N
MASKV(J)=bitxor(MASKV(J),DIRECT(J,L));
VECTOR(J)=MASKV(J)*SCALE;
end
else
%
% GENERATION OF THE POINT FROM "I" AND "N"
%
MASKV=zeros(1,N);
IM=bitxor(I,bitshift(I,-1));
M=0;
while IM~=0 & M<QP,
M=M+1;
if (bitand(IM,1)==1),
for J=1:N
MASKV(J)=bitxor(MASKV(J),DIRECT(J,M));
end
end
IM=bitshift(IM,-1);
end
for J=1:N
VECTOR(J)=MASKV(J)*SCALE;
end
end
%
PRVNUM=I;
PRVDIM=N;
return

14
matlab/cumplot.m Normal file
View File

@ -0,0 +1,14 @@
function h = cumplot(x);
%function h =cumplot(x)
% Copyright (C) 2005 Marco Ratto
n=length(x);
x=[-inf; sort(x); Inf];
y=[0:n n]./n;
h0 = stairs(x,y);
grid on,
if nargout,
h=h0;
end

25
matlab/lptauSEQ.m Normal file
View File

@ -0,0 +1,25 @@
function [lpmat] = lptauSEQ(Nsam,Nvar)
% [lpmat] = lptauSEQ(Nsam,Nvar)
%
% function call LPTAU and generates a sample of dimension Nsam for a
% number of parameters Nvar
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% marco.ratto@jrc.it
%
clear lptau
lpmat = zeros(Nsam, Nvar);
for j=1:Nsam,
lpmat(j,:)=LPTAU(j,Nvar);
end
return

49
matlab/smirnov.m Normal file
View File

@ -0,0 +1,49 @@
function [H,prob,d] = smirnov(x1 , x2 , alpha )
% Smirnov test for 2 distributions
% [H,prob,d] = smirnov(x1 , x2 , alpha )
%
% Copyright (C) 2005 Marco Ratto
%
if nargin<3
alpha = 0.05;
end
% empirical cdfs.
bin = [-inf ; sort([x1;x2]) ; inf];
ncount1 = histc (x1 , bin);
ncount2 = histc (x2 , bin);
cum1 = cumsum(ncount1)./sum(ncount1);
cum2 = cumsum(ncount2)./sum(ncount2);
cum1 = cum1(1:end-1);
cum2 = cum2(1:end-1);
n1 = length(x1);
n2 = length(x2);
n = n1 * n2 /(n1 + n2);
% Compute the d(n1,n2) statistic.
d = max(abs(cum1 - cum2));
%
% Compute P-value check H0 hypothesis
%
lam = max((sqrt(n) + 0.12 + 0.11/sqrt(n)) * d , 0);
j = [1:101]';
prob = 2 * sum((-1).^(j-1).*exp(-2*lam*lam*j.^2));
if prob < 0 , prob = 0; end
if prob > 1 , prob = 1; end
H = (alpha >= prob);

277
matlab/stab_map_.m Normal file
View File

@ -0,0 +1,277 @@
function x0 = stab_map_(Nsam, fload, alpha2, alpha)
%
% function x0 = stab_map_(Nsam, fload, alpha2, alpha)
%
% Mapping of stability regions in the prior ranges applying
% Monte Carlo filtering techniques.
%
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models
% I. Mapping stability, MIMEO, 2005.
%
% INPUTS
% Nsam = MC sample size
% fload = 0 to run new MC; 1 to load prevoiusly generated analysis
% alpha2 = significance level for bivariate sensitivity analysis
% [abs(corrcoef) > alpha2]
% alpha = significance level for univariate sensitivity analysis
% (uses smirnov)
%
% OUTPUT:
% x0: one parameter vector for which the model is stable.
%
% GRAPHS
% 1) Histograms of marginal distributions under the stability regions
% 2) Cumulative distributions of:
% - stable subset (dotted lines)
% - unstable subset (solid lines)
% 3) Bivariate plots of significant correlation patterns
% ( abs(corrcoef) > alpha2) under the stable subset
%
% USES lptauSEQ,
% smirnov
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% marco.ratto@jrc.it
%
% ALL COPIES MUST BE PROVIDED FREE OF CHARGE AND MUST INCLUDE THIS COPYRIGHT
% NOTICE.
%
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
ys_ = oo_.dr.ys;
dr_ = oo_.dr;
fname_ = M_.fname;
nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
if nargin==0,
Nsam=2000; %2^13; %256;
end
if nargin<2,
fload=0;
end
if nargin<4,
alpha=0.002;
end
if nargin<3,
alpha2=0.3;
end
if fload==0 | nargin<2 | isempty(fload),
if estim_params_.np<52,
[lpmat] = lptauSEQ(Nsam,estim_params_.np);
else
%[lpmat] = rand(Nsam,estim_params_.np);
for j=1:estim_params_.np,
lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
end
end
for j=1:estim_params_.np,
if estim_params_.np>30 & estim_params_.np<52
lpmat(:,j)=lpmat(randperm(Nsam),j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
else
lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
end
end
%
h = waitbar(0,'Please wait...');
options_.periods=0;
options_.nomoments=1;
options_.irf=0;
options_.noprint=1;
for j=1:Nsam,
for i=1:estim_params_.np,
evalin('base',[bayestopt_.name{i+nshock}, '= ',sprintf('%0.15e',lpmat(j,i)),';'])
end
%evalin('base','stoch_simul(var_list_);');
stoch_simul([]);
%egg(:,j) = sort(eigenvalues_);
%egg(:,j) = sort(dr_.eigval);
if isfield(dr_,'eigval'),
egg(:,j) = sort(dr_.eigval);
else
egg(:,j)=ones(size(egg,1),1).*1.1;
end
ys_=real(dr_.ys);
yys(:,j) = ys_;
ys_=yys(:,1);
waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
end
close(h)
% map stable samples
ix=[1:Nsam];
for j=1:Nsam,
if abs(egg(dr_.npred+1,j))<=options_.qz_criterium; %1+1.e-5;
ix(j)=0;
elseif abs(egg(dr_.npred,j))>=options_.qz_criterium; %(1-(options_.qz_criterium-1)); %1-1.e-5;
ix(j)=0;
end
end
ix=ix(find(ix)); % stable params
% map unstable samples
ixx=[1:Nsam];
for j=1:Nsam,
%if abs(egg(dr_.npred+1,j))>1+1.e-5 & abs(egg(dr_.npred,j))<1-1.e-5;
if abs(egg(dr_.npred+1,j))>options_.qz_criterium & abs(egg(dr_.npred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
ixx(j)=0;
end
end
ixx=ixx(find(ixx)); % unstable params
save([fname_,'_stab'],'lpmat','ixx','ix','egg','yys')
else
load([fname_,'_stab'])
Nsam = size(lpmat,1);
end
delete([fname_,'_stab_*.*']);
delete([fname_,'_stab_SA_*.*']);
delete([fname_,'_stab_corr_*.*']);
% Blanchard Kahn
for i=1:ceil(estim_params_.np/12),
figure,
for j=1+12*(i-1):min(estim_params_.np,12*i),
subplot(3,4,j-12*(i-1))
hist(lpmat(ix,j),30)
title(bayestopt_.name{j+nshock})
end
saveas(gcf,[fname_,'_stab_',int2str(i)])
end
% TFP STEP & Blanchard;
for i=1:ceil(estim_params_.np/12),
figure,
for j=1+12*(i-1):min(estim_params_.np,12*i),
subplot(3,4,j-12*(i-1))
if ~isempty(ix),
h=cumplot(lpmat(ix,j));
set(h,'color',[1 1 1], 'linestyle',':')
end
hold on,
if ~isempty(ixx),
h=cumplot(lpmat(ixx,j));
set(h,'color',[1 1 1])
end
% if exist('kstest2')==2 & length(ixx)>0 & length(ixx)<Nsam,
% [H,P,KSSTAT] = kstest2(lpmat(ix,j),lpmat(ixx,j));
% title([bayestopt_.name{j+nshock},'. K-S prob ', num2str(P)])
% else
[H,P,KSSTAT] = smirnov(lpmat(ix,j),lpmat(ixx,j));
title([bayestopt_.name{j+nshock},'. K-S prob ', num2str(P)])
% end
end
saveas(gcf,[fname_,'_stab_SA_',int2str(i)])
end
if length(ixx)>0 & length(ixx)<Nsam,
disp(' ')
disp(' ')
disp('Starting bivariate analysis:')
c0=corrcoef(lpmat(ix,:));
c00=tril(c0,-1);
ifig=0;
j2=0;
npar=estim_params_.np;
for j=1:npar,
i2=find(abs(c00(:,j))>alpha2);
if length(i2)>0,
for jx=1:length(i2),
j2=j2+1;
if mod(j2,12)==1,
ifig=ifig+1;
figure('name',['Correlations in the stable sample ', num2str(ifig)]),
end
subplot(3,4,j2-(ifig-1)*12)
% bar(c0(i2,j)),
% set(gca,'xticklabel',bayestopt_.name(i2)),
% set(gca,'xtick',[1:length(i2)])
%plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.w')
%hold on,
plot(lpmat(ix,j),lpmat(ix,i2(jx)),'.')
xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
title(['cc = ',num2str(c0(i2(jx),j))])
if (mod(j2,12)==0) & j2>0,
saveas(gcf,[fname_,'_stab_corr_',int2str(ifig)])
end
end
end
if (j==(npar)) & j2>0,
saveas(gcf,[fname_,'_stab_corr_',int2str(ifig)])
end
end
%close all
end
% % optional map cyclicity of dominant eigenvalues, if
% thex=[];
% for j=1:Nsam,
% %cyc(j)=max(abs(imag(egg(1:34,j))));
% ic = find(imag(egg(1:dr_.npred,j)));
% i=find( abs(egg( ic ,j) )>0.9); %only consider complex dominant eigenvalues
% if ~isempty(i),
% i=i(1:2:end);
% thedum=[];
% for ii=1:length(i),
% idum = ic( i(ii) );
% thedum(ii)=abs(angle(egg(idum,j)));
% end
% [dum, icx]=max(thedum);
% icy(j) = ic( i(icx) );
% thet(j)=max(thedum);
% if thet(j)<0.05 & find(ix==j), % keep stable runs with freq smaller than 0.05
% thex=[thex; j];
% end
% else
% if find(ix==j),
% thex=[thex; j];
% end
% end
% end
% % cyclicity
% for i=1:ceil(estim_params_.np/12),
% figure,
% for j=1+12*(i-1):min(estim_params_.np,12*i),
% subplot(3,4,j-12*(i-1))
% hist(lpmat(thex,j),30)
% title(bayestopt_.name{j+nshock})
% end
% end
%
% % TFP STEP & Blanchard; & cyclicity
% for i=1:ceil(estim_params_.np/12),
% figure,
% for j=1+12*(i-1):min(estim_params_.np,12*i),
% [H,P,KSSTAT] = kstest2(lpmat(1:Nsam,j),lpmat(ixx,j));
% subplot(3,4,j-12*(i-1))
% cdfplot(lpmat(1:Nsam,j))
% hold on,
% cdfplot(lpmat(ixx,j))
% title([bayestopt_.name{j+nshock},'. K-S prob ', num2str(P)])
% end
% end
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
x0 = [x0; lpmat(ix(1),:)'];