-
Notifications
You must be signed in to change notification settings - Fork 0
/
my_World3_sampling_est.m~
1559 lines (1430 loc) · 122 KB
/
my_World3_sampling_est.m~
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
% estimation is carried out via iterative regression wherein each iteration
% obtains parameter estimates given the data at time t and t-1.
% the final parameter estimates come from the mean of the set of parameter
% estimates across all iterations
clear all;
close all;
clc;
addpath(genpath('/home/mrra/CasADi_v3.5.5_linux'));
import casadi.*
%TO DO: fix ML funs so that out of bounds values are more reasonable
% add appropriate variable bounds to the input variables in the tables
skip_ML=false; % choose whether to interpolate and extrapolate the lookup tables provided with the World3 model
%% data
disp('Reading data...')
% load initial states and get variable names
World3_init_guesses = World3_init2();
var_names = fields(World3_init_guesses);
n_var = length(var_names);
% load estimation data
raw_data = readtable('Estimation/my_World3_data.csv');
H = 16; % forecast horizon
% data to use for estimation
est_start = 201; %201 = 1950Q1
est_end = 488-H; %301 = 1975Q1
est_iter_period = 16; % number of timesteps in each iterative regression
minT = max(H,est_iter_period);
dt=0.25; % 0.5 => each time step is 6 months
obs_vars = raw_data(:,4:end).Properties.VariableNames; % skip year, quarter, date columns
% old approach:
% observed variables (MUST MATCH ORDER IN DATA!)
% obs_vars = [...
% "fertility_control_facilities_per_capita";... %1
% "perceived_life_expectancy";... %2
% "delayed_industrial_output_per_capita";... %3
% "desired_completed_family_size";... %4
% "maximum_total_fertility";... %5
% "fecundity_multiplier";... %6
% "population";... %7
% "total_fertility";... %8
% "deaths";... %9
% "Population_0_To_14";... %10
% "Population_15_To_44";... %11
% "Population_45_To_64";... %12
% "Population_65_Plus";...%13
% "food";... %14
% "industrial_output";...%15
% "fraction_of_industrial_output_allocated_to_services_1";... %16
% "Persistent_Pollution";... %17
% "Nonrenewable_Resources";... %18
% "Industrial_Capital";... %19
% "Arable_Land";... %20
% "life_expectancy";... %21
% "inflation";... %22
% "GDP_per_capita";... %23
% ];
% get the variable indices for each observed variable
n_obs = length(obs_vars);
obs_var_inds = nan(n_obs,1);
for v = 1:n_obs
for i = 1:n_var
if strcmp(var_names{i},obs_vars(v)) % column names in data must exactly match variable names, including case
obs_var_inds(v) = i;
break;
end
end
end
% put normalization in a struct and an array for easy look-up
normalization = nan(1,n_var);
for v = 1:n_var
%%% THE PROBLEM HERE when doing so in an automated fashion is that some
%%% quantities may appear together but have different normalizations (ie. x*10^3 - y*10^2)
% norm.(var_names{v}) = 10^floor(log10(max(abs(World3_init_guesses.(var_names{v})),1)));
% normalization(v) = 10^floor(log10(max(abs(World3_init_guesses.(var_names{v})),1)));
if floor(log10(max(abs(World3_init_guesses.(var_names{v})),1)))>6
norm.(var_names{v}) = 10^9;
normalization(v) = 10^9;
else
norm.(var_names{v}) = 1;
normalization(v) = 1;
end
if strcmp(var_names{v},'land_erosion_rate') || strcmp(var_names{v},'land_removal_for_urban_and_industrial_use')
norm.(var_names{v}) = 10^9;
normalization(v) = 10^9;
end
end
% normalize data
data = table2array(raw_data(est_start:4*dt:est_end+H,4:end))./normalization(obs_var_inds); % exclude timestamp column
T=est_end-est_start+1;
infq_data = lin_interp(raw_data.inflation); % inflation data
%% set up shocks and initial conditions
if ~skip_ML
all_shocks = SX.sym('shocks',[n_obs,T]);
shock_counter = 1;
for var = 1:n_var
eval([var_names{var} '_shk = SX.sym(''' var_names{var} '_shk'',[T,1]);']);
% only add shocks/residuals to observed variables
if ~ismember(var,obs_var_inds)
eval([var_names{var} '_shk = zeros(T,1);']);
else
eval(['all_shocks(' num2str(shock_counter) ',:) = ' var_names{var} '_shk;']);
shock_counter = shock_counter + 1;
end
% eval(['all_shocks(' num2str(var) ',:) = ' var_names{var} '_shk;']);
end
end
shocks = cell(T,1);
for t = 1:T
% shocks
shocks{t} = all_shocks(:,t);
end
n_shocks = length(shocks{1});
% % initial state for each variable
% if ~skip_ML
% init_vars = SX.sym('init',[n_var,1]);
% for var = 1:n_var
% eval([var_names{var} '_init = SX.sym(''' var_names{var} '_init '',[1,1]);' ])
% eval(['init_vars(' num2str(var) ') = ' var_names{var} '_init;']); % collect all initial conditions
% % eval([var_names{var} '_init = ' num2str(World3_init_guesses.(var_names{var})) ';']); % assign values to initial conditions
% end
% end
%% Tables:
%%% TO DO: use these tables to create list of upper nad lower bounds on variables
% for example:
% all>0
% 1>frac>0
% potentially_arable_land_total>Potentially_Arable_Land>0
% 1>land_yield_multiplier_from_air_pollution
% 1>land_yield_technology_change_rate_multiplier
% 1>land_life_multiplier_from_land_yield
% 1>marginal_land_yield_multiplier_from_capital
if ~skip_ML
disp('Generating tables...')
names{2}='x'; names{3} = 'y';
x_cas=SX.sym('x',[1,1]);
% Education_Index_LOOKUP=[[0,0];[1000,0.81];[2000,0.88];[3000,0.92];[4000,0.95];[5000,0.98];[6000,0.99];[7000,1]];
% [coeffs{1},fit{1},fun_text{1},modelFun{1}] = get_fit(Education_Index_LOOKUP(:,1),Education_Index_LOOKUP(:,2),false,false,names);
% Education_Index_LOOKUP = @(x) modelFun{1}(coeffs{1},x);
% Education_Index_LOOKUP=Function('f',{x_cas},{Education_Index_LOOKUP(x_cas)});
GDP_per_capita_LOOKUP=[[0,120];[200,600];[400,1200];[600,1800];[800,2500];[1000,3200]];
[coeffs{2},fit{2},fun_text{2},modelFun{2}] = get_fit(GDP_per_capita_LOOKUP(:,1),GDP_per_capita_LOOKUP(:,2),false,false,names);
GDP_per_capita_LOOKUP = @(x) modelFun{2}(coeffs{2},x);
GDP_per_capita_LOOKUP=Function('f',{x_cas},{GDP_per_capita_LOOKUP(x_cas)});
% Life_Expectancy_Index_LOOKUP=[[25,0];[35,0.16];[45,0.33];[55,0.5];[65,0.67];[75,0.84];[85,1]];
% [coeffs{3},fit{3},fun_text{3},modelFun{3}] = get_fit(Life_Expectancy_Index_LOOKUP(:,1),Life_Expectancy_Index_LOOKUP(:,2),false,false,names);
% Life_Expectancy_Index_LOOKUP = @(x) modelFun{3}(coeffs{3},x);
% Life_Expectancy_Index_LOOKUP=Function('f',{x_cas},{Life_Expectancy_Index_LOOKUP(x_cas)});
% development_cost_per_hectare_table=[[0,100000];[0.1,7400];[0.2,5200];[0.3,3500];[0.4,2400];[0.5,1500];[0.6,750];[0.7,300];[0.8,150];[0.9,75];[1,50]];
development_cost_per_hectare_table=[[0,12000];[0.1,7400];[0.2,5200];[0.3,3500];[0.4,2400];[0.5,1500];[0.6,750];[0.7,300];[0.8,150];[0.9,75];[1,50]];
[coeffs{4},fit{4},fun_text{4},modelFun{4}] = get_fit(development_cost_per_hectare_table(:,1),development_cost_per_hectare_table(:,2),false,false,names);
development_cost_per_hectare_table = @(x) modelFun{4}(coeffs{4},x);
development_cost_per_hectare_table=Function('f',{x_cas},{development_cost_per_hectare_table(x_cas)});
% development_cost_per_hectare_table = @(x) interp_f(development_cost_per_hectare_table,x);
% fraction_industrial_output_allocated_to_agriculture_table_1=[[0,0.4];[0.5,0.2];[1,0.1];[1.5,0.025];[2,0];[2.5,0]];
fraction_industrial_output_allocated_to_agriculture_table_1=[[0,0.4];[0.5,0.2];[1,0.1];[1.5,0.025];[2,0.01];[2.5,0.005]];
[coeffs{5},fit{5},fun_text{5},modelFun{5}] = get_fit(fraction_industrial_output_allocated_to_agriculture_table_1(:,1),fraction_industrial_output_allocated_to_agriculture_table_1(:,2),false,false,names);
fraction_industrial_output_allocated_to_agriculture_table_1 = @(x) modelFun{5}(coeffs{5},x);
fraction_industrial_output_allocated_to_agriculture_table_1=Function('f',{x_cas},{fraction_industrial_output_allocated_to_agriculture_table_1(x_cas)});
% fraction_industrial_output_allocated_to_agriculture_table_1 = @(x) interp_f(fraction_industrial_output_allocated_to_agriculture_table_1,x);
% fraction_industrial_output_allocated_to_agriculture_table_2=[[0,0.4];[0.5,0.2];[1,0.1];[1.5,0.025];[2,0];[2.5,0]];
fraction_industrial_output_allocated_to_agriculture_table_2=[[0,0.4];[0.5,0.2];[1,0.1];[1.5,0.025];[2,0.005];[2.5,0]];
[coeffs{6},fit{6},fun_text{6},modelFun{6}] = get_fit(fraction_industrial_output_allocated_to_agriculture_table_2(:,1),fraction_industrial_output_allocated_to_agriculture_table_2(:,2),false,false,names);
fraction_industrial_output_allocated_to_agriculture_table_2 = @(x) modelFun{6}(coeffs{6},x);
fraction_industrial_output_allocated_to_agriculture_table_2=Function('f',{x_cas},{fraction_industrial_output_allocated_to_agriculture_table_2(x_cas)});
% fraction_industrial_output_allocated_to_agriculture_table_2 = @(x) interp_f(fraction_industrial_output_allocated_to_agriculture_table_2,x);
indicated_food_per_capita_table_1=[[0,230];[200,480];[400,690];[600,850];[800,970];[1000,1070];[1200,1150];[1400,1210];[1600,1250]];
[coeffs{7},fit{7},fun_text{7},modelFun{7}] = get_fit(indicated_food_per_capita_table_1(:,1),indicated_food_per_capita_table_1(:,2),false,false,names);
indicated_food_per_capita_table_1 = @(x) modelFun{7}(coeffs{7},x);
indicated_food_per_capita_table_1=Function('f',{x_cas},{indicated_food_per_capita_table_1(x_cas)});
% indicated_food_per_capita_table_1 = @(x) interp_f(indicated_food_per_capita_table_1,x);
indicated_food_per_capita_table_2=[[0,230];[200,480];[400,690];[600,850];[800,970];[1000,1070];[1200,1150];[1400,1210];[1600,1250]];
[coeffs{8},fit{8},fun_text{8},modelFun{8}] = get_fit(indicated_food_per_capita_table_2(:,1),indicated_food_per_capita_table_2(:,2),false,false,names);
indicated_food_per_capita_table_2 = @(x) modelFun{8}(coeffs{8},x);
indicated_food_per_capita_table_2=Function('f',{x_cas},{indicated_food_per_capita_table_2(x_cas)});
% indicated_food_per_capita_table_2 = @(x) interp_f(indicated_food_per_capita_table_2,x);
% land_yield_multiplier_from_capital_table=[[0,1];[40,3];[80,4.5];[120,5];[160,5.3];[200,5.6];[240,5.9];[280,6.1];[320,6.35];[360,6.6];[400,6.9];[440,7.2];[480,7.4];[520,7.6];[560,7.8];[600,8];[640,8.2];[680,8.4];[720,8.6];[760,8.8];[800,9];[840,9.2];[880,9.4];[920,9.6];[960,9.8];[1000,10]];
land_yield_multiplier_from_capital_table=[[0,1];[40,3];[80,3.8];[120,4.4];[160,4.9];[200,5.4];[240,5.7];[280,6.0];[320,6.3];[360,6.6];[400,6.9];[440,7.2];[480,7.4];[520,7.6];[560,7.8];[600,8];[640,8.2];[680,8.4];[720,8.6];[760,8.8];[800,9];[840,9.2];[880,9.4];[920,9.6];[960,9.8];[1000,10]];
% land_yield_multiplier_from_capital_table=[[0,1];[5.33,1.26];[6.576,1.329];[9.96,1.50];[40,3];[80,3.8];[120,4.4];[160,4.9];[200,5.4];[240,5.7];[280,6.0];[320,6.3];[360,6.6];[400,6.9];[440,7.2];[480,7.4];[520,7.6];[560,7.8];[600,8];[640,8.2];[680,8.4];[720,8.6];[760,8.8];[800,9];[840,9.2];[880,9.4];[920,9.6];[960,9.8];[1000,10]];
[coeffs{9},fit{9},fun_text{9},modelFun{9}] = get_fit(land_yield_multiplier_from_capital_table(:,1),land_yield_multiplier_from_capital_table(:,2),false,false,names);
land_yield_multiplier_from_capital_table = @(x) modelFun{9}(coeffs{9},x);
land_yield_multiplier_from_capital_table=Function('f',{x_cas},{land_yield_multiplier_from_capital_table(x_cas)});
% land_yield_multiplier_from_capital_table = @(x) interp_f(land_yield_multiplier_from_capital_table,x);
land_yield_multipler_from_air_pollution_table_1=[[0,1];[10,1];[20,0.7];[30,0.4]];
% land_yield_multipler_from_air_pollution_table_1=[[0,1];[0.5,1];[1,1];[5,1];[10,1];[20,0.7];[30,0.4]];
[coeffs{10},fit{10},fun_text{10},modelFun{10}] = get_fit(land_yield_multipler_from_air_pollution_table_1(:,1),land_yield_multipler_from_air_pollution_table_1(:,2),false,false,names);
land_yield_multipler_from_air_pollution_table_1 = @(x) modelFun{10}(coeffs{10},x);
land_yield_multipler_from_air_pollution_table_1=Function('f',{x_cas},{land_yield_multipler_from_air_pollution_table_1(x_cas)});
% land_yield_multipler_from_air_pollution_table_1 = @(x) interp_f(land_yield_multipler_from_air_pollution_table_1,x);
land_yield_multipler_from_air_pollution_table_2=[[0,1];[10,1];[20,0.98];[30,0.95]];
[coeffs{11},fit{11},fun_text{11},modelFun{11}] = get_fit(land_yield_multipler_from_air_pollution_table_2(:,1),land_yield_multipler_from_air_pollution_table_2(:,2),false,false,names);
land_yield_multipler_from_air_pollution_table_2 = @(x) modelFun{11}(coeffs{11},x);
land_yield_multipler_from_air_pollution_table_2=Function('f',{x_cas},{land_yield_multipler_from_air_pollution_table_2(x_cas)});
% land_yield_multipler_from_air_pollution_table_2 = @(x) interp_f(land_yield_multipler_from_air_pollution_table_2,x);
% land_yield_technology_change_rate_multiplier_table=[[0,0];[1,0]];
land_yield_technology_change_rate_multiplier_table=[[-1,0];[0,0];[1,1];[2,1]];
[coeffs{12},fit{12},fun_text{12},modelFun{12}] = get_fit(land_yield_technology_change_rate_multiplier_table(:,1),land_yield_technology_change_rate_multiplier_table(:,2),false,false,names);
land_yield_technology_change_rate_multiplier_table = @(x) modelFun{12}(coeffs{12},x);
land_yield_technology_change_rate_multiplier_table=Function('f',{x_cas},{land_yield_technology_change_rate_multiplier_table(x_cas)});
% land_yield_technology_change_rate_multiplier_table = @(x) interp_f(land_yield_technology_change_rate_multiplier_table,x);
land_life_multiplier_from_land_yield_table_1=[[0,1.2];[1,1];[2,0.63];[3,0.36];[4,0.16];[5,0.055];[6,0.04];[7,0.025];[8,0.015];[9,0.01]];
% land_life_multiplier_from_land_yield_table_1 = @(x) interp_f(land_life_multiplier_from_land_yield_table_1,x);
[coeffs{13},fit{13},fun_text{13},modelFun{13}] = get_fit(land_life_multiplier_from_land_yield_table_1(:,1),land_life_multiplier_from_land_yield_table_1(:,2),false,false,names);
land_life_multiplier_from_land_yield_table_1 = @(x) modelFun{13}(coeffs{13},x);
land_life_multiplier_from_land_yield_table_1=Function('f',{x_cas},{land_life_multiplier_from_land_yield_table_1(x_cas)});
land_life_multiplier_from_land_yield_table_2=[[0,1.2];[1,1];[2,0.63];[3,0.36];[4,0.29];[5,0.26];[6,0.24];[7,0.22];[8,0.21];[9,0.2]];
% land_life_multiplier_from_land_yield_table_2 = @(x) interp_f(land_life_multiplier_from_land_yield_table_2,x);
[coeffs{14},fit{14},fun_text{14},modelFun{14}] = get_fit(land_life_multiplier_from_land_yield_table_2(:,1),land_life_multiplier_from_land_yield_table_2(:,2),false,false,names);
land_life_multiplier_from_land_yield_table_2 = @(x) modelFun{14}(coeffs{14},x);
land_life_multiplier_from_land_yield_table_2=Function('f',{x_cas},{land_life_multiplier_from_land_yield_table_2(x_cas)});
urban_and_industrial_land_required_per_capita_table=[[0,0.005];[200,0.008];[400,0.015];[600,0.025];[800,0.04];[1000,0.055];[1200,0.07];[1400,0.08];[1600,0.09]];
% urban_and_industrial_land_required_per_capita_table = @(x) interp_f(urban_and_industrial_land_required_per_capita_table,x);
[coeffs{15},fit{15},fun_text{15},modelFun{15}] = get_fit(urban_and_industrial_land_required_per_capita_table(:,1),urban_and_industrial_land_required_per_capita_table(:,2),false,false,names);
urban_and_industrial_land_required_per_capita_table = @(x) modelFun{15}(coeffs{15},x);
urban_and_industrial_land_required_per_capita_table=Function('f',{x_cas},{urban_and_industrial_land_required_per_capita_table(x_cas)});
land_fertility_degredation_rate_table=[[0,0];[10,0.1];[20,0.3];[30,0.5]];
% land_fertility_degredation_rate_table = @(x) interp_f(land_fertility_degredation_rate_table,x);
[coeffs{16},fit{16},fun_text{16},modelFun{16}] = get_fit(land_fertility_degredation_rate_table(:,1),land_fertility_degredation_rate_table(:,2),false,false,names);
coeffs{16}(3)=0; %fixed point at 0, and always positive for positive x
land_fertility_degredation_rate_table = @(x) modelFun{16}(coeffs{16},x);
land_fertility_degredation_rate_table=Function('f',{x_cas},{land_fertility_degredation_rate_table(x_cas)});
land_fertility_regeneration_time_table=[[0,20];[0.02,13];[0.04,8];[0.06,4];[0.08,2];[0.1,2]];
% land_fertility_regeneration_time_table = @(x) interp_f(land_fertility_regeneration_time_table,x);
[coeffs{17},fit{17},fun_text{17},modelFun{17}] = get_fit(land_fertility_regeneration_time_table(:,1),land_fertility_regeneration_time_table(:,2),false,false,names);
land_fertility_regeneration_time_table = @(x) modelFun{17}(coeffs{17},x);
land_fertility_regeneration_time_table=Function('f',{x_cas},{land_fertility_regeneration_time_table(x_cas)});
fraction_of_agricultural_inputs_for_land_maintenance_table=[[0,0];[1,0.04];[2,0.07];[3,0.09];[4,0.1]];
% fraction_of_agricultural_inputs_for_land_maintenance_table = @(x) interp_f(fraction_of_agricultural_inputs_for_land_maintenance_table,x);
[coeffs{18},fit{18},fun_text{18},modelFun{18}] = get_fit(fraction_of_agricultural_inputs_for_land_maintenance_table(:,1),fraction_of_agricultural_inputs_for_land_maintenance_table(:,2),false,false,names);
fraction_of_agricultural_inputs_for_land_maintenance_table = @(x) modelFun{18}(coeffs{18},x);
fraction_of_agricultural_inputs_for_land_maintenance_table=Function('f',{x_cas},{fraction_of_agricultural_inputs_for_land_maintenance_table(x_cas)});
fraction_of_agricultural_inputs_allocated_to_land_dev_table=[[0,0];[0.25,0.05];[0.5,0.15];[0.75,0.3];[1,0.5];[1.25,0.7];[1.5,0.85];[1.75,0.95];[2,1]];
% fraction_of_agricultural_inputs_allocated_to_land_dev_table = @(x) interp_f(fraction_of_agricultural_inputs_allocated_to_land_dev_table,x);
[coeffs{19},fit{19},fun_text{19},modelFun{19}] = get_fit(fraction_of_agricultural_inputs_allocated_to_land_dev_table(:,1),fraction_of_agricultural_inputs_allocated_to_land_dev_table(:,2),false,false,names);
fraction_of_agricultural_inputs_allocated_to_land_dev_table = @(x) modelFun{19}(coeffs{19},x);
fraction_of_agricultural_inputs_allocated_to_land_dev_table=Function('f',{x_cas},{fraction_of_agricultural_inputs_allocated_to_land_dev_table(x_cas)});
marginal_land_yield_multiplier_from_capital_table=[[0,0.075];[40,0.03];[80,0.015];[120,0.011];[160,0.009];[200,0.008];[240,0.007];[280,0.006];[320,0.005];[360,0.005];[400,0.005];[440,0.005];[480,0.005];[520,0.005];[560,0.005];[600,0.005]];
% marginal_land_yield_multiplier_from_capital_table = @(x) interp_f(marginal_land_yield_multiplier_from_capital_table,x);
[coeffs{20},fit{20},fun_text{20},modelFun{20}] = get_fit(marginal_land_yield_multiplier_from_capital_table(:,1),marginal_land_yield_multiplier_from_capital_table(:,2),false,false,names);
marginal_land_yield_multiplier_from_capital_table = @(x) modelFun{20}(coeffs{20},x);
marginal_land_yield_multiplier_from_capital_table=Function('f',{x_cas},{marginal_land_yield_multiplier_from_capital_table(x_cas)});
industrial_capital_output_ratio_multiplier_from_resource_table=[[0,3.75];[0.1,3.6];[0.2,3.47];[0.3,3.36];[0.4,3.25];[0.5,3.16];[0.6,3.1];[0.7,3.06];[0.8,3.02];[0.9,3.01];[1,3]];
% industrial_capital_output_ratio_multiplier_from_resource_table = @(x) interp_f(industrial_capital_output_ratio_multiplier_from_resource_table,x);
[coeffs{21},fit{21},fun_text{21},modelFun{21}] = get_fit(industrial_capital_output_ratio_multiplier_from_resource_table(:,1),industrial_capital_output_ratio_multiplier_from_resource_table(:,2),false,false,names);
industrial_capital_output_ratio_multiplier_from_resource_table = @(x) modelFun{21}(coeffs{21},x);
industrial_capital_output_ratio_multiplier_from_resource_table=Function('f',{x_cas},{industrial_capital_output_ratio_multiplier_from_resource_table(x_cas)});
frac_of_industrial_output_allocated_to_consumption_var_table=[[0,0.3];[0.2,0.32];[0.4,0.34];[0.6,0.36];[0.8,0.38];[1,0.43];[1.2,0.73];[1.4,0.77];[1.6,0.81];[1.8,0.82];[2,0.83]];
% frac_of_industrial_output_allocated_to_consumption_var_table = @(x) interp_f(frac_of_industrial_output_allocated_to_consumption_var_table,x);
[coeffs{22},fit{22},fun_text{22},modelFun{22}] = get_fit(frac_of_industrial_output_allocated_to_consumption_var_table(:,1),frac_of_industrial_output_allocated_to_consumption_var_table(:,2),false,false,names);
frac_of_industrial_output_allocated_to_consumption_var_table = @(x) modelFun{22}(coeffs{22},x);
frac_of_industrial_output_allocated_to_consumption_var_table=Function('f',{x_cas},{frac_of_industrial_output_allocated_to_consumption_var_table(x_cas)});
industrial_capital_output_ratio_multiplier_from_pollution_table=[[0,1.25];[0.1,1.2];[0.2,1.15];[0.3,1.11];[0.4,1.08];[0.5,1.05];[0.6,1.03];[0.7,1.02];[0.8,1.01];[0.9,1];[1,1]];
% industrial_capital_output_ratio_multiplier_from_pollution_table = @(x) interp_f(industrial_capital_output_ratio_multiplier_from_pollution_table,x);
[coeffs{23},fit{23},fun_text{23},modelFun{23}] = get_fit(industrial_capital_output_ratio_multiplier_from_pollution_table(:,1),industrial_capital_output_ratio_multiplier_from_pollution_table(:,2),false,false,names);
industrial_capital_output_ratio_multiplier_from_pollution_table = @(x) modelFun{23}(coeffs{23},x);
industrial_capital_output_ratio_multiplier_from_pollution_table=Function('f',{x_cas},{industrial_capital_output_ratio_multiplier_from_pollution_table(x_cas)});
industrial_capital_output_ratio_multiplier_table=[[1,1];[1.2,1.05];[1.4,1.12];[1.6,1.25];[1.8,1.35];[2,1.5]];
% industrial_capital_output_ratio_multiplier_table = @(x) interp_f(industrial_capital_output_ratio_multiplier_table,x);
[coeffs{24},fit{24},fun_text{24},modelFun{24}] = get_fit(industrial_capital_output_ratio_multiplier_table(:,1),industrial_capital_output_ratio_multiplier_table(:,2),false,false,names);
industrial_capital_output_ratio_multiplier_table = @(x) modelFun{24}(coeffs{24},x);
industrial_capital_output_ratio_multiplier_table=Function('f',{x_cas},{industrial_capital_output_ratio_multiplier_table(x_cas)});
capacity_utilization_fraction_table=[[1,1];[3,0.9];[5,0.7];[7,0.3];[9,0.1];[11,0.1]];
% capacity_utilization_fraction_table = @(x) interp_f(capacity_utilization_fraction_table,x);
[coeffs{25},fit{25},fun_text{25},modelFun{25}] = get_fit(capacity_utilization_fraction_table(:,1),capacity_utilization_fraction_table(:,2),false,false,names);
capacity_utilization_fraction_table = @(x) modelFun{25}(coeffs{25},x);
capacity_utilization_fraction_table=Function('f',{x_cas},{capacity_utilization_fraction_table(x_cas)});
jobs_per_hectare_table=[[2,2];[6,0.5];[10,0.4];[14,0.3];[18,0.27];[22,0.24];[26,0.2];[30,0.2]];
% jobs_per_hectare_table=[[0,5];[2,2];[6,0.5];[10,0.4];[14,0.3];[18,0.27];[22,0.24];[30,0.2];[50,0.15]];
% jobs_per_hectare_table = @(x) interp_f(jobs_per_hectare_table,x);
[coeffs{26},fit{26},fun_text{26},modelFun{26}] = get_fit(jobs_per_hectare_table(:,1),jobs_per_hectare_table(:,2),false,false,names);
jobs_per_hectare_table = @(x) modelFun{26}(coeffs{26},x);
jobs_per_hectare_table=Function('f',{x_cas},{jobs_per_hectare_table(x_cas)});
jobs_per_industrial_capital_unit_table=[[50,0.37];[200,0.18];[350,0.12];[500,0.09];[650,0.07];[800,0.06]];
% jobs_per_industrial_capital_unit_table = @(x) interp_f(jobs_per_industrial_capital_unit_table,x);
[coeffs{27},fit{27},fun_text{27},modelFun{27}] = get_fit(jobs_per_industrial_capital_unit_table(:,1),jobs_per_industrial_capital_unit_table(:,2),false,false,names);
jobs_per_industrial_capital_unit_table = @(x) modelFun{27}(coeffs{27},x);
jobs_per_industrial_capital_unit_table=Function('f',{x_cas},{jobs_per_industrial_capital_unit_table(x_cas)});
jobs_per_service_capital_unit_table=[[50,1.1];[200,0.6];[350,0.35];[500,0.2];[650,0.15];[800,0.15];[950,0.15];[1100,0.15];[1250,0.15];[1400,0.15];[1550,0.15]];
% jobs_per_service_capital_unit_table = @(x) interp_f(jobs_per_service_capital_unit_table,x);
[coeffs{28},fit{28},fun_text{28},modelFun{28}] = get_fit(jobs_per_service_capital_unit_table(:,1),jobs_per_service_capital_unit_table(:,2),false,false,names);
jobs_per_service_capital_unit_table = @(x) modelFun{28}(coeffs{28},x);
jobs_per_service_capital_unit_table=Function('f',{x_cas},{jobs_per_service_capital_unit_table(x_cas)});
fraction_of_industrial_output_allocated_to_services_table_1=[[0,0.3];[0.5,0.2];[1,0.1];[1.5,0.05];[2,0]];
% fraction_of_industrial_output_allocated_to_services_table_1 = @(x) interp_f(fraction_of_industrial_output_allocated_to_services_table_1,x);
[coeffs{29},fit{29},fun_text{29},modelFun{29}] = get_fit(fraction_of_industrial_output_allocated_to_services_table_1(:,1),fraction_of_industrial_output_allocated_to_services_table_1(:,2),false,false,names);
fraction_of_industrial_output_allocated_to_services_table_1 = @(x) modelFun{29}(coeffs{29},x);
fraction_of_industrial_output_allocated_to_services_table_1=Function('f',{x_cas},{fraction_of_industrial_output_allocated_to_services_table_1(x_cas)});
fraction_of_industrial_output_allocated_to_services_table_2=[[0,0.3];[0.5,0.2];[1,0.1];[1.5,0.05];[2,0]];
% fraction_of_industrial_output_allocated_to_services_table_2 = @(x) interp_f(fraction_of_industrial_output_allocated_to_services_table_2,x);
[coeffs{30},fit{30},fun_text{30},modelFun{30}] = get_fit(fraction_of_industrial_output_allocated_to_services_table_2(:,1),fraction_of_industrial_output_allocated_to_services_table_2(:,2),false,false,names);
fraction_of_industrial_output_allocated_to_services_table_2 = @(x) modelFun{30}(coeffs{30},x);
fraction_of_industrial_output_allocated_to_services_table_2=Function('f',{x_cas},{fraction_of_industrial_output_allocated_to_services_table_2(x_cas)});
indicated_services_output_per_capita_table_1=[[0,40];[200,300];[400,640];[600,1000];[800,1220];[1000,1450];[1200,1650];[1400,1800];[1600,2000]];
% indicated_services_output_per_capita_table_1=[[0,40];[41.5,94];[200,300];[400,640];[600,1000];[800,1220];[1000,1450];[1200,1650];[1400,1800];[1600,2000]];
% indicated_services_output_per_capita_table_1 = @(x) interp_f(indicated_services_output_per_capita_table_1,x);
[coeffs{31},fit{31},fun_text{31},modelFun{31}] = get_fit(indicated_services_output_per_capita_table_1(:,1),indicated_services_output_per_capita_table_1(:,2),false,false,names);
indicated_services_output_per_capita_table_1 = @(x) modelFun{31}(coeffs{31},x);
indicated_services_output_per_capita_table_1=Function('f',{x_cas},{indicated_services_output_per_capita_table_1(x_cas)});
indicated_services_output_per_capita_table_2=[[0,40];[200,300];[400,640];[600,1000];[800,1220];[1000,1450];[1200,1650];[1400,1800];[1600,2000]];
% indicated_services_output_per_capita_table_2 = @(x) interp_f(indicated_services_output_per_capita_table_2,x);
[coeffs{32},fit{32},fun_text{32},modelFun{32}] = get_fit(indicated_services_output_per_capita_table_2(:,1),indicated_services_output_per_capita_table_2(:,2),false,false,names);
indicated_services_output_per_capita_table_2 = @(x) modelFun{32}(coeffs{32},x);
indicated_services_output_per_capita_table_2=Function('f',{x_cas},{indicated_services_output_per_capita_table_2(x_cas)});
assimilation_half_life_mult_table=[[1,1];[251,11];[501,21];[751,31];[1001,41]];
% assimilation_half_life_mult_table=[[0,1];[0.5,1];[1,1];[2.29,1.05];[3.53,1.10];[7.60,1.26];[9.439,1.333];[11,1.4];[12,1.44];[13,1.48];[14,1.52];[15,1.56];[251,11];[501,21];[751,31];[1001,41]];
% assimilation_half_life_mult_table = @(x) interp_f(assimilation_half_life_mult_table,x);
[coeffs{33},fit{33},fun_text{33},modelFun{33}] = get_fit(assimilation_half_life_mult_table(:,1),assimilation_half_life_mult_table(:,2),false,false,names);
assimilation_half_life_mult_table = @(x) modelFun{33}(coeffs{33},x);
assimilation_half_life_mult_table=Function('f',{x_cas},{assimilation_half_life_mult_table(x_cas)});
% persistent_pollution_technology_change_mult_table=[[-1,0];[0,0]];
persistent_pollution_technology_change_mult_table=[[-1,0];[0,0];[1,0]];
% persistent_pollution_technology_change_mult_table = @(x) interp_f(persistent_pollution_technology_change_mult_table,x);
[coeffs{34},fit{34},fun_text{34},modelFun{34}] = get_fit(persistent_pollution_technology_change_mult_table(:,1),persistent_pollution_technology_change_mult_table(:,2),false,false,names);
persistent_pollution_technology_change_mult_table = @(x) modelFun{34}(coeffs{34},x);
persistent_pollution_technology_change_mult_table=Function('f',{x_cas},{persistent_pollution_technology_change_mult_table(x_cas)});
mortality_45_to_64_table=[[20,0.0562];[30,0.0373];[40,0.0252];[50,0.0171];[60,0.0118];[70,0.0083];[80,0.006]];
% mortality_45_to_64_table = @(x) interp_f(mortality_45_to_64_table,x);
[coeffs{35},fit{35},fun_text{35},modelFun{35}] = get_fit(mortality_45_to_64_table(:,1),mortality_45_to_64_table(:,2),false,false,names);
mortality_45_to_64_table = @(x) modelFun{35}(coeffs{35},x);
mortality_45_to_64_table=Function('f',{x_cas},{mortality_45_to_64_table(x_cas)});
mortality_65_plus_table=[[20,0.13];[30,0.11];[40,0.09];[50,0.07];[60,0.06];[70,0.05];[80,0.04]];
% mortality_65_plus_table = @(x) interp_f(mortality_65_plus_table,x);
[coeffs{37},fit{37},fun_text{37},modelFun{37}] = get_fit(mortality_65_plus_table(:,1),mortality_65_plus_table(:,2),false,false,names);
mortality_65_plus_table = @(x) modelFun{37}(coeffs{37},x);
mortality_65_plus_table=Function('f',{x_cas},{mortality_65_plus_table(x_cas)});
mortality_0_to_14_table=[[20,0.0567];[30,0.0366];[40,0.0243];[50,0.0155];[60,0.0082];[70,0.0023];[80,0.001]];
% mortality_0_to_14_table = @(x) interp_f(mortality_0_to_14_table,x);
[coeffs{39},fit{39},fun_text{39},modelFun{39}] = get_fit(mortality_0_to_14_table(:,1),mortality_0_to_14_table(:,2),false,false,names);
mortality_0_to_14_table = @(x) modelFun{39}(coeffs{39},x);
mortality_0_to_14_table=Function('f',{x_cas},{mortality_0_to_14_table(x_cas)});
mortality_15_to_44_table=[[20,0.0266];[30,0.0171];[40,0.011];[50,0.0065];[60,0.004];[70,0.0016];[80,0.0008]];
% mortality_15_to_44_table = @(x) interp_f(mortality_15_to_44_table,x);
[coeffs{41},fit{41},fun_text{41},modelFun{41}] = get_fit(mortality_15_to_44_table(:,1),mortality_15_to_44_table(:,2),false,false,names);
mortality_15_to_44_table = @(x) modelFun{41}(coeffs{41},x);
mortality_15_to_44_table=Function('f',{x_cas},{mortality_15_to_44_table(x_cas)});
completed_multiplier_from_perceived_lifetime_table=[[0,3];[10,2.1];[20,1.6];[30,1.4];[40,1.3];[50,1.2];[60,1.1];[70,1.05];[80,1]];
% completed_multiplier_from_perceived_lifetime_table = @(x) interp_f(completed_multiplier_from_perceived_lifetime_table,x);
[coeffs{43},fit{43},fun_text{43},modelFun{43}] = get_fit(completed_multiplier_from_perceived_lifetime_table(:,1),completed_multiplier_from_perceived_lifetime_table(:,2),false,false,names);
completed_multiplier_from_perceived_lifetime_table = @(x) modelFun{43}(coeffs{43},x);
completed_multiplier_from_perceived_lifetime_table=Function('f',{x_cas},{completed_multiplier_from_perceived_lifetime_table(x_cas)});
family_response_to_social_norm_table=[[-0.2,0.5];[-0.1,0.6];[0,0.7];[0.1,0.85];[0.2,1]];
% family_response_to_social_norm_table = @(x) interp_f(family_response_to_social_norm_table,x);
[coeffs{44},fit{44},fun_text{44},modelFun{44}] = get_fit(family_response_to_social_norm_table(:,1),family_response_to_social_norm_table(:,2),false,false,names);
family_response_to_social_norm_table = @(x) modelFun{44}(coeffs{44},x);
family_response_to_social_norm_table=Function('f',{x_cas},{family_response_to_social_norm_table(x_cas)});
fecundity_multiplier_table=[[0,0];[10,0.2];[20,0.4];[30,0.6];[40,0.7];[50,0.75];[60,0.79];[70,0.84];[80,0.87]];
% fecundity_multiplier_table = @(x) interp_f(fecundity_multiplier_table,x);
[coeffs{46},fit{46},fun_text{46},modelFun{46}] = get_fit(fecundity_multiplier_table(:,1),fecundity_multiplier_table(:,2),false,false,names);
fecundity_multiplier_table = @(x) modelFun{46}(coeffs{46},x);
fecundity_multiplier_table=Function('f',{x_cas},{fecundity_multiplier_table(x_cas)});
fertility_control_effectiveness_table=[[0,0.75];[0.5,0.85];[1,0.9];[1.5,0.95];[2,0.98];[2.5,0.99];[3,1]];
% fertility_control_effectiveness_table=[[0,0.75];[0.5,0.85];[1,0.9];[1.5,0.93];[2,0.95];[2.5,0.96];[3,0.97];[4,0.98];[5,0.99];[6,0.993];[7,0.993];[8,0.995];[9,0.997];[10,0.999]];
% fertility_control_effectiveness_table = @(x) interp_f(fertility_control_effectiveness_table,x);
[coeffs{47},fit{47},fun_text{47},modelFun{47}] = get_fit(fertility_control_effectiveness_table(:,1),fertility_control_effectiveness_table(:,2),false,false,names);
fertility_control_effectiveness_table = @(x) modelFun{47}(coeffs{47},x);
fertility_control_effectiveness_table=Function('f',{x_cas},{fertility_control_effectiveness_table(x_cas)});
fraction_services_allocated_to_fertility_control_table=[[-1,0];[0,0];[2,0.005];[4,0.015];[6,0.025];[8,0.03];[10,0.035];[20,0.04];[30,0.045]];
% fraction_services_allocated_to_fertility_control_table = @(x) interp_f(fraction_services_allocated_to_fertility_control_table,x);
[coeffs{49},fit{49},fun_text{49},modelFun{49}] = get_fit(fraction_services_allocated_to_fertility_control_table(:,1),fraction_services_allocated_to_fertility_control_table(:,2),false,false,names);
fraction_services_allocated_to_fertility_control_table = @(x) modelFun{49}(coeffs{49},x);
fraction_services_allocated_to_fertility_control_table=Function('f',{x_cas},{fraction_services_allocated_to_fertility_control_table(x_cas)});
% social_family_size_normal_table=[[0,1.25];[200,0.94];[400,0.715];[600,0.59];[800,0.5]];
social_family_size_normal_table=[[0,1.25];[200,1];[400,0.9];[600,0.8];[800,0.75]];
% social_family_size_normal_table = @(x) interp_f(social_family_size_normal_table,x);
[coeffs{50},fit{50},fun_text{50},modelFun{50}] = get_fit(social_family_size_normal_table(:,1),social_family_size_normal_table(:,2),false,false,names);
social_family_size_normal_table = @(x) modelFun{50}(coeffs{50},x);
social_family_size_normal_table=Function('f',{x_cas},{social_family_size_normal_table(x_cas)});
crowding_multiplier_from_industry_table=[[0,0.5];[200,0.05];[400,-0.1];[600,-0.08];[800,-0.02];[1000,0.05];[1200,0.1];[1400,0.15];[1600,0.2]];
% crowding_multiplier_from_industry_table = @(x) interp_f(crowding_multiplier_from_industry_table,x);
[coeffs{51},fit{51},fun_text{51},modelFun{51}] = get_fit(crowding_multiplier_from_industry_table(:,1),crowding_multiplier_from_industry_table(:,2),false,false,names);
crowding_multiplier_from_industry_table = @(x) modelFun{51}(coeffs{51},x);
crowding_multiplier_from_industry_table=Function('f',{x_cas},{crowding_multiplier_from_industry_table(x_cas)});
fraction_of_population_urban_table=[[0,0];[2e+009,0.2];[4e+009,0.4];[6e+009,0.5];[8e+009,0.58];[1e+010,0.65];[1.2e+010,0.72];[1.4e+010,0.78];[1.6e+010,0.8]];
% fraction_of_population_urban_table = @(x) interp_f(fraction_of_population_urban_table,x);
[coeffs{52},fit{52},fun_text{52},modelFun{52}] = get_fit(fraction_of_population_urban_table(:,1),fraction_of_population_urban_table(:,2),false,false,names);
fraction_of_population_urban_table = @(x) modelFun{52}(coeffs{52},x);
fraction_of_population_urban_table=Function('f',{x_cas},{fraction_of_population_urban_table(x_cas)});
% health_services_per_capita_table=[[0,0];[250,20];[500,50];[750,95];[1000,140];[1250,175];[1500,200];[1750,220];[2000,230]];
health_services_per_capita_table=[[0,0];[90,7.2];[250,20];[500,50];[750,95];[1000,140];[1250,175];[1500,200];[1750,220];[2000,230]];
% health_services_per_capita_table = @(x) interp_f(health_services_per_capita_table,x);
[coeffs{53},fit{53},fun_text{53},modelFun{53}] = get_fit(health_services_per_capita_table(:,1),health_services_per_capita_table(:,2),false,false,names);
health_services_per_capita_table = @(x) modelFun{53}(coeffs{53},x);
health_services_per_capita_table=Function('f',{x_cas},{health_services_per_capita_table(x_cas)});
% lifetime_multiplier_from_food_table=[[0,0];[1,1];[2,1.43];[3,1.5];[4,1.5];[5,1.5]];
lifetime_multiplier_from_food_table=[[0,0];[1,1];[2,1.2];[3,1.3];[4,1.35];[5,1.4]];
% lifetime_multiplier_from_food_table = @(x) interp_f(lifetime_multiplier_from_food_table,x);
[coeffs{55},fit{55},fun_text{55},modelFun{55}] = get_fit(lifetime_multiplier_from_food_table(:,1),lifetime_multiplier_from_food_table(:,2),false,false,names);
lifetime_multiplier_from_food_table = @(x) modelFun{55}(coeffs{55},x);
lifetime_multiplier_from_food_table=Function('f',{x_cas},{lifetime_multiplier_from_food_table(x_cas)});
lifetime_multiplier_from_health_services_1_table=[[0,1];[20,1.1];[40,1.4];[60,1.6];[80,1.7];[100,1.8]];
% lifetime_multiplier_from_health_services_1_table = @(x) interp_f(lifetime_multiplier_from_health_services_1_table,x);
[coeffs{57},fit{57},fun_text{57},modelFun{57}] = get_fit(lifetime_multiplier_from_health_services_1_table(:,1),lifetime_multiplier_from_health_services_1_table(:,2),false,false,names);
lifetime_multiplier_from_health_services_1_table = @(x) modelFun{57}(coeffs{57},x);
lifetime_multiplier_from_health_services_1_table=Function('f',{x_cas},{lifetime_multiplier_from_health_services_1_table(x_cas)});
% lifetime_multiplier_from_health_services_2_table=[[0,1];[20,1.5];[40,1.9];[60,2];[80,2];[100,2]];
lifetime_multiplier_from_health_services_2_table=[[0,1];[20,1.4];[40,1.6];[60,1.8];[80,1.95];[100,2]];
% lifetime_multiplier_from_health_services_2_table = @(x) interp_f(lifetime_multiplier_from_health_services_2_table,x);
[coeffs{58},fit{58},fun_text{58},modelFun{58}] = get_fit(lifetime_multiplier_from_health_services_2_table(:,1),lifetime_multiplier_from_health_services_2_table(:,2),false,false,names);
lifetime_multiplier_from_health_services_2_table = @(x) modelFun{58}(coeffs{58},x);
lifetime_multiplier_from_health_services_2_table=Function('f',{x_cas},{lifetime_multiplier_from_health_services_2_table(x_cas)});
lifetime_multiplier_from_persistent_pollution_table=[[0,1];[10,0.99];[20,0.97];[30,0.95];[40,0.9];[50,0.85];[60,0.75];[70,0.65];[80,0.55];[90,0.4];[100,0.2]];
% lifetime_multiplier_from_persistent_pollution_table=[[0,1];[2,0.999];[4,0.998];[6,0.997];[8,0.996];[10,0.99];[20,0.97];[30,0.95];[40,0.9];[50,0.85];[60,0.75];[70,0.65];[80,0.55];[90,0.4];[100,0.2];[200,0.01];[300,0.01];[400,0.001]];
% lifetime_multiplier_from_persistent_pollution_table = @(x) interp_f(lifetime_multiplier_from_persistent_pollution_table,x);
[coeffs{59},fit{59},fun_text{59},modelFun{59}] = get_fit(lifetime_multiplier_from_persistent_pollution_table(:,1),lifetime_multiplier_from_persistent_pollution_table(:,2),false,false,names);
lifetime_multiplier_from_persistent_pollution_table = @(x) modelFun{59}(coeffs{59},x);
lifetime_multiplier_from_persistent_pollution_table=Function('f',{x_cas},{lifetime_multiplier_from_persistent_pollution_table(x_cas)});
fraction_of_capital_allocated_to_obtaining_resources_1_table=[[0,1];[0.1,0.9];[0.2,0.7];[0.3,0.5];[0.4,0.2];[0.5,0.1];[0.6,0.05];[0.7,0.05];[0.8,0.05];[0.9,0.05];[1,0.05]];
% fraction_of_capital_allocated_to_obtaining_resources_1_table = @(x) interp_f(fraction_of_capital_allocated_to_obtaining_resources_1_table,x);
[coeffs{60},fit{60},fun_text{60},modelFun{60}] = get_fit(fraction_of_capital_allocated_to_obtaining_resources_1_table(:,1),fraction_of_capital_allocated_to_obtaining_resources_1_table(:,2),false,false,names);
fraction_of_capital_allocated_to_obtaining_resources_1_table = @(x) modelFun{60}(coeffs{60},x);
fraction_of_capital_allocated_to_obtaining_resources_1_table=Function('f',{x_cas},{fraction_of_capital_allocated_to_obtaining_resources_1_table(x_cas)});
fraction_of_capital_allocated_to_obtaining_resources_2_table=[[0,1];[0.1,0.2];[0.2,0.1];[0.3,0.05];[0.4,0.05];[0.5,0.05];[0.6,0.05];[0.7,0.05];[0.8,0.05];[0.9,0.05];[1,0.05]];
% fraction_of_capital_allocated_to_obtaining_resources_2_table = @(x) interp_f(fraction_of_capital_allocated_to_obtaining_resources_2_table,x);
[coeffs{61},fit{61},fun_text{61},modelFun{61}] = get_fit(fraction_of_capital_allocated_to_obtaining_resources_2_table(:,1),fraction_of_capital_allocated_to_obtaining_resources_2_table(:,2),false,false,names);
fraction_of_capital_allocated_to_obtaining_resources_2_table = @(x) modelFun{61}(coeffs{61},x);
fraction_of_capital_allocated_to_obtaining_resources_2_table=Function('f',{x_cas},{fraction_of_capital_allocated_to_obtaining_resources_2_table(x_cas)});
% resource_technology_change_mult_table=[[-1,0];[0,0]];
resource_technology_change_mult_table=[[-1,0];[0,0];[1,0]];
% resource_technology_change_mult_table = @(x) interp_f(resource_technology_change_mult_table,x);
[coeffs{62},fit{62},fun_text{62},modelFun{62}] = get_fit(resource_technology_change_mult_table(:,1),resource_technology_change_mult_table(:,2),false,false,names);
resource_technology_change_mult_table = @(x) modelFun{62}(coeffs{62},x);
resource_technology_change_mult_table=Function('f',{x_cas},{resource_technology_change_mult_table(x_cas)});
% per_capita_resource_use_mult_table=[[0,0];[200,0.85];[400,2.6];[600,3.4];[800,3.8];[1000,4.1];[1200,4.4];[1400,4.7];[1600,5]];
per_capita_resource_use_mult_table=[[0,0];[200,0.85];[400,2.6];[600,4.4];[800,5.4];[1000,6.2];[1200,6.8];[1400,7.0];[1600,7.0]];
% per_capita_resource_use_mult_table=[[0,0];[20,0.085];[40,0.17];[100,0.43];[200,0.85];[400,2.6];[600,4.4];[800,5.4];[1000,6.2];[1200,6.8];[1400,7.0];[1600,7.0]];
% per_capita_resource_use_mult_table = @(x) interp_f(per_capita_resource_use_mult_table,x);
[coeffs{63},fit{63},fun_text{63},modelFun{63}] = get_fit(per_capita_resource_use_mult_table(:,1),per_capita_resource_use_mult_table(:,2),false,false,names);
per_capita_resource_use_mult_table = @(x) modelFun{63}(coeffs{63},x);
per_capita_resource_use_mult_table=Function('f',{x_cas},{per_capita_resource_use_mult_table(x_cas)});
end
%% Calibrated Params:
disp('Generating parameters...')
persistent_pollution_transmission_delay=20;
% arable_land_erosion_factor=0.92;
% w3_real_exhange_rate=7.18783;
LIFE_TIME_MULTIPLIER_FROM_SERVICES=1;%0.557;
% investment_into_services=1;%2.07;
GDP_pc_unit=1;
unit_agricultural_input=1;
unit_population=1;
ha_per_Gha=1e+009;
% ha_per_unit_of_pollution=4;
one_year=1;
% Ref_Hi_GDP=9508;
% Ref_Lo_GDP=24;
% Total_Land=1.91;
% land_fraction_harvested=0.7;%0.68;
potentially_arable_land_total=3.2e+009/norm.('Potentially_Arable_Land');
% processing_loss=0.1;%0.35;
% desired_food_ratio=2;
IND_OUT_IN_1970=7.9e+011/norm.('industrial_output');
average_life_of_agricultural_inputs_1=2;
average_life_of_agricultural_inputs_2=2;
% land_yield_factor_1=1;%1.634;
% air_pollution_policy_implementation_time=4000;
% average_life_of_land_normal=6000/norm.('average_life_of_land');%1000;
% land_life_policy_implementation_time=4000;
urban_and_industrial_land_development_time=10;
% inherent_land_fertility=600;
food_shortage_perception_delay=2;
subsistence_food_per_capita=230/norm.('food_per_capita');
% social_discount=0.07;
% industrial_output_per_capita_desired=400;
average_life_of_industrial_capital_1=14/norm.('average_life_of_industrial_capital');
average_life_of_industrial_capital_2=14/norm.('average_life_of_industrial_capital');
fraction_of_industrial_output_allocated_to_consumption_const_1=0.43;
fraction_of_industrial_output_allocated_to_consumption_const_2=0.43;
% industrial_capital_output_ratio_1=3;%3.64;
% industrial_equilibrium_time=4000;
labor_force_participation_fraction=0.75;
labor_utilization_fraction_delay_time=2;
average_life_of_service_capital_1=24.07/norm.('average_life_of_service_capital');
average_life_of_service_capital_2=20/norm.('average_life_of_service_capital');
% service_capital_output_ratio_1=1;%0.579;
service_capital_output_ratio_2=1;
% FINAL_TIME=2100;
% INITIAL_TIME=1995;
% SAVEPER=1;
% POLICY_YEAR=2101;
% TIME_STEP=0.0078125;
agricultural_material_toxicity_index=1;
assimilation_half_life_in_1970=1.5/norm.('assimilation_half_life');%0.95;
% desired_persistent_pollution_index=1.2;
% fraction_of_agricultural_inputs_from_persistent_materials=0.001;
% fraction_of_resources_from_persistent_materials=0.02;
% industrial_material_toxicity_index=10;
% industrial_material_emissions_factor=0.1;
% persistent_pollution_generation_factor_1=1;%0;
persistent_pollution_in_1970=1.36e+008/norm.('Persistent_Pollution');%2.5e7;
% desired_completed_family_size_normal=4;%3.412;
income_expectation_averaging_time=3;%3.572;
lifetime_perception_delay=20;%23.54;
% maximum_total_fertility_normal=12;
reproductive_lifetime=30;
social_adjustment_delay=20;%26;
% fertility_control_effectiveness_time=4000;
% population_equilibrium_time=4000;
% zero_population_growth_time=4000;
% THOUSAND=1000;
health_services_impact_delay=20;%24.977;
% life_expectancy_normal=28;
desired_resource_use_rate=4.8e+009/norm.('resource_usage_rate');
% resource_use_factor_1=1;%0;
% frac_of_ind_capital_allocated_to_obtaining_res_switch_time=1990;
technology_development_delay=20;
% PRICE_OF_FOOD=0.22;
%% params to estimate
% maximum_total_fertility_normal = SX.sym('maximum_total_fertility_normal',1,1);
% desired_completed_family_size_normal = SX.sym('desired_completed_family_size_normal',1,1);
% land_fraction_harvested = SX.sym('land_fraction_harvested',1,1);
% desired_food_ratio = SX.sym('desired_food_ratio',1,1);
% inherent_land_fertility = SX.sym('inherent_land_fertility',1,1);
% inv_industrial_output_per_capita_desired = SX.sym('industrial_output_per_capita_desired',1,1);
% inv_desired_persistent_pollution_index = SX.sym('desired_persistent_pollution_index',1,1);
% inv_social_discount = SX.sym('social_discount',1,1);
% industrial_capital_output_ratio_1 = SX.sym('industrial_capital_output_ratio_1',1,1);
% investment_into_services = SX.sym('investment_into_services',1,1);
% processing_loss = SX.sym('processing_loss',1,1);
% average_life_of_land_normal = SX.sym('average_life_of_land_normal',1,1);
% fraction_of_agricultural_inputs_from_persistent_materials = SX.sym('fraction_of_agricultural_inputs_from_persistent_materials',1,1);
% fraction_of_resources_from_persistent_materials = SX.sym('fraction_of_resources_from_persistent_materials',1,1);
% industrial_material_toxicity_index = SX.sym('industrial_material_toxicity_index',1,1);
% industrial_material_emissions_factor = SX.sym('industrial_material_emissions_factor',1,1);
% persistent_pollution_generation_rate_factor = SX.sym('persistent_pollution_generation_rate_factor',1,1);
% life_expectancy_normal = SX.sym('life_expectancy_normal',1,1);
% resource_use_factor_1 = SX.sym('resource_use_factor_1',1,1);
% persistent_pollution_generation_factor_1 = SX.sym('persistent_pollution_generation_factor_1',1,1);
% land_yield_factor_1 = SX.sym('land_yield_factor_1',1,1);
% service_capital_output_ratio_1 = SX.sym('service_capital_output_ratio_1',1,1);
%
% fcfpc_param = SX.sym('fcfpc_param',1,1);
% le_param1 = SX.sym('le_param1',1,1);
% le_param2 = SX.sym('le_param2',1,1);
% diopc_param1 = SX.sym('diopc_param1',1,1);
% diopc_param2 = SX.sym('diopc_param2',1,1);
% ppar_param1 = SX.sym('ppar_param1',1,1);
% ppar_param2 = SX.sym('ppar_param2',1,1);
% ppar_param3 = SX.sym('ppar_param3',1,1);
%
% inflation_param1 = SX.sym('inflation_param1',1,1);
% inflation_param2 = SX.sym('inflation_param2',1,1);
% inflation_param3 = SX.sym('inflation_param3',1,1);
% inflation_param4 = SX.sym('inflation_param4',1,1);
% inflation_param5 = SX.sym('inflation_param5',1,1);
% inflation_param6 = SX.sym('inflation_param6',1,1);
% inflation_init_1 = SX.sym('inflation_init_1',1,1);
% inflation_init_2 = SX.sym('inflation_init_2',1,1);
% inflation_init_3 = SX.sym('inflation_init_3',1,1);
% inflation_init_4 = SX.sym('inflation_init_4',1,1);
param_names = ["maximum_total_fertility_normal";...
"desired_completed_family_size_normal";...
"land_fraction_harvested";...
"desired_food_ratio";...
"inherent_land_fertility";...
"inv_industrial_output_per_capita_desired";...
"inv_desired_persistent_pollution_index";...
"inv_social_discount";...
"industrial_capital_output_ratio_1";...
"investment_into_services";...
"processing_loss";...
"average_life_of_land_normal";...
"fraction_of_agricultural_inputs_from_persistent_materials";...
"fraction_of_resources_from_persistent_materials";...
"industrial_material_toxicity_index";...
"industrial_material_emissions_factor";...
"persistent_pollution_generation_rate_factor";...
"life_expectancy_normal";...
"fcfpc_param";...
"le_param1";...
"le_param2";...
"diopc_param1";...
"diopc_param2";...
"resource_use_factor_1";...
"ppar_param1";...
"ppar_param2";...
"ppar_param3";...
"persistent_pollution_generation_factor_1";...
"land_yield_factor_1";...
"service_capital_output_ratio_1";...
"inflation_param1";...
"inflation_param2";...
"inflation_param3";...
"inflation_param4";...
"inflation_param5";...
"inflation_param6";...
% "inflation_init_1";...
% "inflation_init_2";...
% "inflation_init_3";...
% "inflation_init_4";...
];
n_params = length(param_names);
params = SX.sym('params',n_params,1);
for p = 1:n_params
eval([char(param_names(p)) ' = SX.sym(''' char(param_names(p)) ''',1,1);']);
eval(['params(p) = ' char(param_names(p)) ';']);
end
%%% TO DO: get proper initial guesses for:
% fcfpc_param
% inherent_land_fertility
%%% TO DO: put this in a separate file and then load it
param_settings.lower_bound.maximum_total_fertility_normal = 8;%2;
param_settings.init_param_guess.maximum_total_fertility_normal = 12;
param_settings.upper_bound.maximum_total_fertility_normal = 16;%20;
param_settings.lower_bound.desired_completed_family_size_normal = 3;%1;
param_settings.init_param_guess.desired_completed_family_size_normal = 3.3;%4;
param_settings.upper_bound.desired_completed_family_size_normal = 5;%8;
param_settings.lower_bound.land_fraction_harvested = 0.4;
param_settings.init_param_guess.land_fraction_harvested = 0.7;
param_settings.upper_bound.land_fraction_harvested = 1;
param_settings.lower_bound.desired_food_ratio = 1;
param_settings.init_param_guess.desired_food_ratio = 2;
param_settings.upper_bound.desired_food_ratio = 3;%5;
param_settings.lower_bound.inherent_land_fertility = 5;
param_settings.init_param_guess.inherent_land_fertility = 6;
param_settings.upper_bound.inherent_land_fertility = 10;
param_settings.lower_bound.inv_industrial_output_per_capita_desired = 0.0001;
param_settings.init_param_guess.inv_industrial_output_per_capita_desired = 1/400;
param_settings.upper_bound.inv_industrial_output_per_capita_desired = 2/400;%0.01;
param_settings.lower_bound.inv_desired_persistent_pollution_index = 1/1.2;%0.5;
param_settings.init_param_guess.inv_desired_persistent_pollution_index = 1/1.2;
param_settings.upper_bound.inv_desired_persistent_pollution_index = 1/1.2;%1;
param_settings.lower_bound.inv_social_discount = 14;%1;
param_settings.init_param_guess.inv_social_discount = 1/0.07;
param_settings.upper_bound.inv_social_discount = 14.5;%50;
param_settings.lower_bound.industrial_capital_output_ratio_1 = 2.95;%1;
param_settings.init_param_guess.industrial_capital_output_ratio_1 = 3;
param_settings.upper_bound.industrial_capital_output_ratio_1 = 3.05;%10;
param_settings.lower_bound.investment_into_services = 0.1;
param_settings.init_param_guess.investment_into_services = 1;
param_settings.upper_bound.investment_into_services = 3;
param_settings.lower_bound.processing_loss = 0;
param_settings.init_param_guess.processing_loss = 0.1;
param_settings.upper_bound.processing_loss = 0.2;
param_settings.lower_bound.average_life_of_land_normal = 5.5;%1;
param_settings.init_param_guess.average_life_of_land_normal = 6;
param_settings.upper_bound.average_life_of_land_normal = 6.5;%10;
param_settings.lower_bound.fraction_of_agricultural_inputs_from_persistent_materials = 0.001;%0.0001;
param_settings.init_param_guess.fraction_of_agricultural_inputs_from_persistent_materials = 0.001;
param_settings.upper_bound.fraction_of_agricultural_inputs_from_persistent_materials = 0.001;%0.1;
param_settings.lower_bound.fraction_of_resources_from_persistent_materials = 0.02;%0.0001;
param_settings.init_param_guess.fraction_of_resources_from_persistent_materials = 0.02;
param_settings.upper_bound.fraction_of_resources_from_persistent_materials = 0.02;%0.1;
param_settings.lower_bound.industrial_material_toxicity_index = 9;%0.1;
param_settings.init_param_guess.industrial_material_toxicity_index = 10;
param_settings.upper_bound.industrial_material_toxicity_index = 11;%100;
param_settings.lower_bound.industrial_material_emissions_factor = 0.1;%0.0001;
param_settings.init_param_guess.industrial_material_emissions_factor = 0.1;
param_settings.upper_bound.industrial_material_emissions_factor = 0.1;%1;
param_settings.lower_bound.persistent_pollution_generation_rate_factor = (0.5*3/persistent_pollution_transmission_delay)^3;%10^-5;
param_settings.init_param_guess.persistent_pollution_generation_rate_factor = (0.5*3/persistent_pollution_transmission_delay)^3;
param_settings.upper_bound.persistent_pollution_generation_rate_factor = (0.5*3/persistent_pollution_transmission_delay)^3;%10^-3;
param_settings.lower_bound.life_expectancy_normal = 20;
param_settings.init_param_guess.life_expectancy_normal = 28;
param_settings.upper_bound.life_expectancy_normal = 50;
param_settings.lower_bound.fcfpc_param = 0;
param_settings.init_param_guess.fcfpc_param = 0.002;
param_settings.upper_bound.fcfpc_param = 5;
param_settings.lower_bound.le_param1 = 0;
param_settings.init_param_guess.le_param1 = 0.02;
param_settings.upper_bound.le_param1 = 5;
param_settings.lower_bound.le_param2 = 0;
param_settings.init_param_guess.le_param2 = 1;
param_settings.upper_bound.le_param2 = 5;
param_settings.lower_bound.diopc_param1 = 0;
param_settings.init_param_guess.diopc_param1 = 0.915;
param_settings.upper_bound.diopc_param1 = 5;
param_settings.lower_bound.diopc_param2 = 0;
param_settings.init_param_guess.diopc_param2 = 1.05;
param_settings.upper_bound.diopc_param2 = 5;
param_settings.lower_bound.resource_use_factor_1 = 0.9;%0;
param_settings.init_param_guess.resource_use_factor_1=1.1;%1;%0.61;
param_settings.upper_bound.resource_use_factor_1 = 1.2;%5;
param_settings.lower_bound.ppar_param1 = (1 + 1/(-5 + 1/(5 + 1/(-8 - 1/3))));%0.7415;%-5;
param_settings.init_param_guess.ppar_param1 = (1 + 1/(-5 + 1/(5 + 1/(-8 - 1/3))));
param_settings.upper_bound.ppar_param1 = (1 + 1/(-5 + 1/(5 + 1/(-8 - 1/3))));%0.8;%5;
param_settings.lower_bound.ppar_param2 = (-3 + 1/(2 + 1/(3 + 1/(4 + 1/(5 + 1/10)))));%-2.5769;%-5;
param_settings.init_param_guess.ppar_param2 = (-3 + 1/(2 + 1/(3 + 1/(4 + 1/(5 + 1/10)))));
param_settings.upper_bound.ppar_param2 = (-3 + 1/(2 + 1/(3 + 1/(4 + 1/(5 + 1/10)))));%-2.5569;%5;
param_settings.lower_bound.ppar_param3 = (0.5*3*(1+3*12)/20);%2.7650;%-5;
param_settings.init_param_guess.ppar_param3 = (0.5*3*(1+3*12)/20);
param_settings.upper_bound.ppar_param3 = (0.5*3*(1+3*12)/20);%2.7850;%5;
param_settings.lower_bound.persistent_pollution_generation_factor_1 = 1.017;%0.95;%0;
param_settings.init_param_guess.persistent_pollution_generation_factor_1 = 1.017;%1;
param_settings.upper_bound.persistent_pollution_generation_factor_1 = 1.017;%1.05;%5;
param_settings.lower_bound.land_yield_factor_1 = 1;%0;
param_settings.init_param_guess.land_yield_factor_1 = 1;%1.86;
param_settings.upper_bound.land_yield_factor_1 = 1.1;%5;
param_settings.lower_bound.service_capital_output_ratio_1 = 0.9;%0;
param_settings.init_param_guess.service_capital_output_ratio_1 = 1;%0.56;
param_settings.upper_bound.service_capital_output_ratio_1 = 1.1;%5;
param_settings.lower_bound.inflation_param1 = -10;
param_settings.init_param_guess.inflation_param1 = 1;
param_settings.upper_bound.inflation_param1 = 10;
param_settings.lower_bound.inflation_param2 = -10;
param_settings.init_param_guess.inflation_param2 = 0;
param_settings.upper_bound.inflation_param2 = 10;
param_settings.lower_bound.inflation_param3 = -10;
param_settings.init_param_guess.inflation_param3 = 0;
param_settings.upper_bound.inflation_param3 = 10;
param_settings.lower_bound.inflation_param4 = -10;
param_settings.init_param_guess.inflation_param4 = 0;
param_settings.upper_bound.inflation_param4 = 10;
param_settings.lower_bound.inflation_param5 = -1;
param_settings.init_param_guess.inflation_param5 = 0;
param_settings.upper_bound.inflation_param5 = 10;
param_settings.lower_bound.inflation_param6 = -1;
param_settings.init_param_guess.inflation_param6 = 0;
param_settings.upper_bound.inflation_param6 = 10;
% param_settings.lower_bound.inflation_init_1 = -10;
% param_settings.init_param_guess.inflation_init_1 = 5.7;
% param_settings.upper_bound.inflation_init_1 = 20;
%
% param_settings.lower_bound.inflation_init_2 = -10;
% param_settings.init_param_guess.inflation_init_2 = 5.7;
% param_settings.upper_bound.inflation_init_2 = 20;
%
% param_settings.lower_bound.inflation_init_3 = -10;
% param_settings.init_param_guess.inflation_init_3 = 5.7;
% param_settings.upper_bound.inflation_init_3 = 20;
%
% param_settings.lower_bound.inflation_init_4 = -10;
% param_settings.init_param_guess.inflation_init_4 = 5.7;
% param_settings.upper_bound.inflation_init_4 = 20;
% initial parameter guess
try
load('param_mean.mat');
init_param_guess = param_means;
disp("loading initial guess from mat file...")
catch
init_param_guess = nan(n_params,1);
for p = 1:n_params
init_param_guess(p) = param_settings.init_param_guess.(param_names(p));
end
init_param_guess = init_param_guess;
end
% parameter lower bounds
param_lb = nan(n_params,1);
for p = 1:n_params
param_lb(p) = param_settings.lower_bound.(param_names(p));
end
param_lb = param_lb;
% parameter upper bounds
param_ub = nan(n_params,1);
for p = 1:n_params
param_ub(p) = param_settings.upper_bound.(param_names(p));
end
param_ub = param_ub;
%% initial conditions are for 1950. If dt = 0.25, a shift=-40 implies a shift in 1940, 100=1975
shift_all = T+1;
shift_fraction_of_industrial_output_allocated_to_agriculture=shift_all;
shift_indicated_food_per_capita=shift_all;
shift_average_life_agricultural_inputs=shift_all;
% shift_resource_technology_change_rate=max(0,(301-est_start)/(dt*4));%shift_all;
shift_fraction_of_industrial_capital_for_obtaining_resources=shift_all;
% shift_resource_use_factor=max(0,(341-est_start)/(dt*4));%shift_all;
% shift_lifetime_multiplier_from_health_services=max(0,(161-est_start)/(dt*4));
shift_fertility_control_effectiveness=T*10;%shift_all;
shift_births=T*10;%shift_all;
shift_desired_completed_family_size=T*10;%shift_all;
shift_persistent_pollution_technology_change_rate=shift_all;
% shift_persistent_pollution_generation_factor=max(0,(321-est_start)/(dt*4));%shift_all;
shift_average_life_of_service_capital=0;%shift_all;
shift_fraction_of_industrial_output_allocated_to_services=shift_all;
shift_indicated_services_output_per_capita=shift_all;
shift_service_capital_output_ratio=shift_all;
shift_average_life_of_industrial_capital=shift_all;
shift_fraction_of_industrial_output_allocated_to_consumption=shift_all;
shift_industrial_capital_output_ratio=shift_all;
shift_fraction_of_industrial_output_for_consumption_const=shift_all;
shift_land_yield_multiplier_from_technology=shift_all;
shift_land_yield_multiplier_from_air_pollution=shift_all;
shift_land_yield_technology_change_rate=shift_all;
shift_land_life_multiplier_from_land_yield=shift_all;
shift_resource_technology_change_rate=SX.sym('shift_resource_technology_change_rate',[1,1]);
shift_resource_use_factor=SX.sym('shift_resource_use_factor',[1,1]);
shift_lifetime_multiplier_from_health_services=SX.sym('shift_lifetime_multiplier_from_health_services',[1,1]);
shift_persistent_pollution_generation_factor=SX.sym('shift_persistent_pollution_generation_factor',[1,1]);
shifts = [shift_resource_technology_change_rate;...
shift_resource_use_factor;...
shift_lifetime_multiplier_from_health_services;...
shift_persistent_pollution_generation_factor];
shift_times = [(301-est_start)/(dt*4);...
(341-est_start)/(dt*4);...
(161-est_start)/(dt*4);...
(321-est_start)/(dt*4)];
%% load a simulation
% this is needed to initialize unobserved state variables
load("world3_sim.mat"); % load the simulation trajectory (state variable values are needed for this estimation method)
if isstruct(sim_traj)
sim_traj = struct2array(sim_traj)';
end
if size(sim_traj,1)~=n_var
error('a simulated trajectory is not available for all variables. Check that your simulation (such as World3f) has all of the required variables and that its output has been saved to world3_sim.mat')
end
if dt~=0.5 % assume that sim_traj has been generated with dt=0.5, since the equations seem to work best with that step length for some reason
sim_traj_temp = nan(n_var,T+H);
sim_traj_temp(:,1:max(1,1/(2*dt)):T+H) = sim_traj(:,est_start-201+1:max(1,dt*2):(est_start-201+T+H)*max(1,dt*2)/max(1,1/(2*dt)));
sim_traj = lin_interp(sim_traj_temp')';
end
%% Dynamics:
disp('Generating equations...')
eq = cell(minT,1);
ineq2eq = cell(minT,1);
% non_nan_data = fill_nan_vals(data);
non_nan_data = lin_interp(data);
% initialize equations
for var = 1:n_var
eval([var_names{var} ' = cell(minT,1);']);
end
init_vars = SX.sym('init',[n_var,1]);
for var = 1:n_var
eval([var_names{var} '_init = SX.sym(''' var_names{var} '_init'',[1,1]);' ])
eval(['init_vars(' num2str(var) ') = ' var_names{var} '_init;']); % collect all initial conditions
end
inflation_init_1 = SX.sym(')
%%% TO DO: add forestry and desertification, perceived vs actual NRR
for t = 1:minT
% Population
if t>1
% eq{t} [
Population_0_To_14{t}=Population_0_To_14_shk(t)+(Population_0_To_14{t-1}+dt*((births{t-1}-deaths_0_to_14{t-1}-maturation_14_to_15{t-1})));...
Population_15_To_44{t}=Population_15_To_44_shk(t)+(Population_15_To_44{t-1}+dt*((maturation_14_to_15{t-1}-deaths_15_to_44{t-1}-maturation_44_to_45{t-1})));...
Population_45_To_64{t}=Population_45_To_64_shk(t)+(Population_45_To_64{t-1}+dt*((maturation_44_to_45{t-1}-deaths_45_to_64{t-1}-maturation_64_to_65{t-1})));...
Population_65_Plus{t}=Population_65_Plus_shk(t)+(Population_65_Plus{t-1}+dt*((maturation_64_to_65{t-1}-deaths_65_plus{t-1})));...
% ];
else
% eq{t} [
% Population_0_To_14{t}=Population_0_To_14_shk(t)+(Population_0_To_14_init+dt*((births_init-deaths_0_to_14_init-maturation_14_to_15_init)));...
% Population_15_To_44{t}=Population_15_To_44_shk(t)+(Population_15_To_44_init+dt*((maturation_14_to_15_init-deaths_15_to_44_init-maturation_44_to_45_init)));...
% Population_45_To_64{t}=Population_45_To_64_shk(t)+(Population_45_To_64_init+dt*((maturation_44_to_45_init-deaths_45_to_64_init-maturation_64_to_65_init)));...
% Population_65_Plus{t}=Population_65_Plus_shk(t)+(Population_65_Plus_init+dt*((maturation_64_to_65_init-deaths_65_plus_init)));...
Population_0_To_14{t}=Population_0_To_14_shk(t)+Population_0_To_14_init;...
Population_15_To_44{t}=Population_15_To_44_shk(t)+Population_15_To_44_init;...
Population_45_To_64{t}=Population_45_To_64_shk(t)+Population_45_To_64_init;...
Population_65_Plus{t}=Population_65_Plus_shk(t)+Population_65_Plus_init;...
% ];
end
% eq{t} [eq{t};...
population{t}=population_shk(t)+(Population_0_To_14{t}+Population_15_To_44{t}+Population_45_To_64{t}+Population_65_Plus{t});...
% ];
% Resources 1
if t>1
% eq{t} [eq{t};...
Nonrenewable_Resources{t}=Nonrenewable_Resources_shk(t)+(Nonrenewable_Resources{t-1}+dt*((-resource_usage_rate{t-1})));...
Resource_Conservation_Technology{t}=Resource_Conservation_Technology_shk(t)+(Resource_Conservation_Technology{t-1}+dt*(resource_technology_change_rate{t-1}));
% resource_use_fact_2{t}=resource_use_fact_2_shk(t)+((resource_use_fact_2{t-1}+dt*sum(vertcat(Resource_Conservation_Technology{max(1,t-technology_development_delay):t-1}))/min(t-1,technology_development_delay))/2);...
resource_use_fact_2{t}=0.7;
% ];
else
% eq{t} [eq{t};...
% Nonrenewable_Resources{t}=Nonrenewable_Resources_shk(t)+(Nonrenewable_Resources_init+dt*((-resource_usage_rate_init)));...
% Resource_Conservation_Technology{t}=Resource_Conservation_Technology_shk(t)+(Resource_Conservation_Technology_init+dt*(resource_technology_change_rate_init));
% resource_use_fact_2{t}=resource_use_fact_2_shk(t)+((resource_use_fact_2_init)/2);...
Nonrenewable_Resources{t}=Nonrenewable_Resources_shk(t)+Nonrenewable_Resources_init;...
Resource_Conservation_Technology{t}=Resource_Conservation_Technology_shk(t)+Resource_Conservation_Technology_init;
resource_use_fact_2{t}=0.7;%resource_use_fact_2_shk(t)+resource_use_fact_2_init;...
% ];
end
% eq{t} [eq{t};...
resource_use_factor{t}=resource_use_factor_shk(t)+(resource_use_fact_2{t}/(1+exp(-40*(t-shift_resource_use_factor)))-resource_use_factor_1/(1+exp(-40*(t-shift_resource_use_factor)))+resource_use_factor_1);...
% ];
fraction_of_resources_remaining{t}=fraction_of_resources_remaining_shk(t)+(Nonrenewable_Resources{t}/Nonrenewable_Resources_init);
% Agriculture 1
if t>1
% eq{t} [eq{t};...
land_yield_factor_2{t}=land_yield_factor_2_shk(t)+((land_yield_factor_2{t-1}+dt*sum(vertcat(Land_Yield_Technology{max(1,t-technology_development_delay):t-1}))/min(t-1,technology_development_delay))/2);...
% ];
else
% eq{t} [eq{t};...
% land_yield_factor_2{t}=land_yield_factor_2_shk(t)+((land_yield_factor_2_init)/2);...
land_yield_factor_2{t}=land_yield_factor_2_shk(t)+land_yield_factor_2_init;...
% ];
end
% eq{t} [eq{t};...
land_yield_multiplier_from_technology{t}=land_yield_multiplier_from_technology_shk(t)+(land_yield_factor_2{t}/(1+exp(-40*(t-shift_land_yield_multiplier_from_technology)))-land_yield_factor_1/(1+exp(-40*(t-shift_land_yield_multiplier_from_technology)))+land_yield_factor_1);...
% ];
% Pollution 1
if t>1
% eq{t} [eq{t};...
Persistent_Pollution_Technology{t}=Persistent_Pollution_Technology_shk(t)+(Persistent_Pollution_Technology{t-1}+dt*(persistent_pollution_technology_change_rate{t-1}));...
% persistent_pollution_generation_factor_2{t}=persistent_pollution_generation_factor_2_shk(t)+((persistent_pollution_generation_factor_2{t-1}+dt*sum(vertcat(Persistent_Pollution_Technology{max(1,t-technology_development_delay):t-1}))/min(t-1,technology_development_delay))/2);...
persistent_pollution_generation_factor_2{t}=0.7;
Persistent_Pollution{t}=Persistent_Pollution_shk(t)+(Persistent_Pollution{t-1}+dt*((persistent_pollution_appearance_rate{t-1}-persistent_pollution_assimilation_rate{t-1})));...
% ];
else
% eq{t} [eq{t};...
% Persistent_Pollution_Technology{t}=Persistent_Pollution_Technology_shk(t)+(Persistent_Pollution_Technology_init+dt*(persistent_pollution_technology_change_rate_init));...
% persistent_pollution_generation_factor_2{t}=persistent_pollution_generation_factor_2_shk(t)+((persistent_pollution_generation_factor_2_init)/2);...
% Persistent_Pollution{t}=Persistent_Pollution_shk(t)+(Persistent_Pollution_init+dt*((persistent_pollution_appearance_rate_init-persistent_pollution_assimilation_rate_init)));...
Persistent_Pollution_Technology{t}=Persistent_Pollution_Technology_shk(t)+Persistent_Pollution_Technology_init;...
persistent_pollution_generation_factor_2{t}=0.7;%persistent_pollution_generation_factor_2_shk(t)+persistent_pollution_generation_factor_2_init;...
Persistent_Pollution{t}=Persistent_Pollution_shk(t)+Persistent_Pollution_init;...
% ];
end
% eq{t} [eq{t};...
persistent_pollution_generation_factor{t}=persistent_pollution_generation_factor_shk(t)+(persistent_pollution_generation_factor_2{t}/(1+exp(-40*(t-shift_persistent_pollution_generation_factor)))-persistent_pollution_generation_factor_1/(1+exp(-40*(t-shift_persistent_pollution_generation_factor)))+persistent_pollution_generation_factor_1);...
% ];
persistent_pollution_index{t}=persistent_pollution_index_shk(t)+(Persistent_Pollution{t}/persistent_pollution_in_1970);
% Industry 1
if t>1
% eq{t} [eq{t};...
Industrial_Capital{t}=Industrial_Capital_shk(t)+(Industrial_Capital{t-1}+dt*((industrial_capital_investment{t-1}-industrial_capital_depreciation{t-1})));...
Delayed_Labor_Utilization_Fraction{t}=Delayed_Labor_Utilization_Fraction_shk(t)+(Delayed_Labor_Utilization_Fraction{t-1}+dt*(labor_utilization_fraction{t-1}-Delayed_Labor_Utilization_Fraction{t-1})/labor_utilization_fraction_delay_time);...
Service_Capital{t}=Service_Capital_shk(t)+(Service_Capital{t-1}+dt*((service_capital_investment{t-1}-service_capital_depreciation{t-1})));...
% ];
else
% eq{t} [eq{t};...
% Industrial_Capital{t}=Industrial_Capital_shk(t)+(Industrial_Capital_init+dt*((industrial_capital_investment_init-industrial_capital_depreciation_init)));...
% Delayed_Labor_Utilization_Fraction{t}=Delayed_Labor_Utilization_Fraction_shk(t)+(Delayed_Labor_Utilization_Fraction_init+dt*(labor_utilization_fraction_init-Delayed_Labor_Utilization_Fraction_init)/labor_utilization_fraction_delay_time);...