REM REM STDP_ex2_new REM Peter Lohmander 171112 CLS OPEN "PLEX2_new.txt" FOR OUTPUT AS #2 PRINT "" PRINT " CCF harvest optimization via stochastic dynamic programming" PRINT " Peter Lohmander Example 2 2017-11-12" PRINT "" PRINT #2, "" PRINT #2, " CCF harvest optimization via stochastic dynamic programming" PRINT #2, " Peter Lohmander Example 2 2017-11-12" PRINT #2, "" DIM f(100, 10, 15), u(100, 10, 15), DPProb(10), xlevel(15), deltaxlevel(15) DIM harvset(15, 15) REM Delta Price Probability distribution DPProb(0) = .00 DPProb(1) = .03 DPProb(2) = .07 DPProb(3) = .10 DPProb(4) = .18 DPProb(5) = .24 DPProb(6) = .18 DPProb(7) = .10 DPProb(8) = .07 DPProb(9) = .03 DPProb(10) = .00 tmax = 5 deltat = 10 r = 0.03 p0 = 10 p1 = 0.3 p2 = 10 c = 50 Landv = 50 REM Transition matrix determination REM http://www.lohmander.com/Moscow_2010/Lohmander_Zazykina_Moscow_2010.ppt REM dx/dt = alfa*x-beta*xx = (1/20)x-(1/8000)xx a = 0.05 alfa = 1 / 20 beta = 1 / 8000 REM K = h = alfa / beta K = alfa / beta x0 = 30 xlevel(0) = x0 deltaxlevel(0) = 0 FOR level = 1 TO 15 xlevel(level) = 1 / (1 / K + (1 / xlevel(level - 1) - 1 / K) * EXP(-a * deltat)) deltaxlevel(level) = xlevel(level) - xlevel(level - 1) NEXT level REM GOTO 1 PRINT " level xlevel deltaxlevel" FOR level = 0 TO 15 PRINT USING "##########"; level; xlevel(level); deltaxlevel(level) NEXT level INPUT zzz 1 REM REM Harvest set determination FOR stocklevel = 0 TO 15 FOR harvlevel = 0 TO 15 harvset(stocklevel, harvlevel) = 0 NEXT harvlevel NEXT stocklevel FOR stocklevel = 1 TO 15 FOR harvlevel = 0 TO stocklevel harvset(stocklevel, harvlevel) = xlevel(stocklevel) - xlevel(stocklevel - harvlevel) NEXT harvlevel NEXT stocklevel GOTO 2 FOR stocklevel = 1 TO 15 PRINT "stocklevel = "; stocklevel; "stock = "; xlevel(stocklevel) PRINT "Harvest alternatives = "; FOR harvlevel = 0 TO stocklevel PRINT USING "####"; harvset(stocklevel, harvlevel); NEXT harvlevel PRINT "" PRINT "" NEXT stocklevel INPUT zzz 2 REM FOR t = 0 TO 100 FOR pricelevel = 0 TO 9 FOR stocklevel = 0 TO 15 f(t, pricelevel, stocklevel) = 0 u(t, pricelevel, stocklevel) = 0 NEXT stocklevel NEXT pricelevel NEXT t REM boundary conditions at tmax t = tmax * deltat FOR pricelevel = 0 TO 9 netp = p0 + p1 * t + p2 * (pricelevel - 5) FOR stocklevel = 0 TO 15 vol = xlevel(stocklevel) f(tmax, pricelevel, stocklevel) = EXP(-r * t) * (netp * vol + Landv) u(tmax, pricelevel, stocklevel) = vol IF f(tmax, pricelevel, stocklevel) < 0 THEN u(tmax, pricelevel, stocklevel) = 0 IF f(tmax, pricelevel, stocklevel) < 0 THEN f(tmax, pricelevel, stocklevel) = 0 NEXT stocklevel NEXT pricelevel REM PRINT "" REM PRINT " f(tmax,pricelevel (=column), stocklevel (=row))" REM PRINT "" REM FOR stocklevel = 1 TO 15 REM FOR pricelevel = 1 TO 9 REM PRINT USING "######"; f(tmax, pricelevel, stocklevel); REM NEXT pricelevel REM PRINT "" REM NEXT stocklevel REM INPUT zzz REM ************************************************************************************************************************* REM Stochastic dynamic programming via backward recursion FOR tindex = tmax - 1 TO 0 STEP -1 t = tindex * deltat disc = EXP(-r * t) FOR stocklevel = 1 TO 15 FOR pricelevel = 1 TO 9 netp = p0 + p1 * t + p2 * (pricelevel - 5) fopt = -999999 uopt = 0 FOR harvlevel = 0 TO stocklevel stocklevel2 = stocklevel + 1 - harvlevel IF stocklevel2 > 15 THEN GOTO 1001 harv = harvset(stocklevel, harvlevel) fnow = EXP(-r * t) * (netp * harv - c) ffut = 0 FOR j = 1 TO 9 ffut = ffut + DPProb(j) * f(tindex + 1, j, stocklevel + 1 - harvlevel) NEXT j fev = fnow + ffut IF fev > fopt THEN uopt = harv IF fev > fopt THEN fopt = fev 1001 REM NEXT harvlevel f(tindex, pricelevel, stocklevel) = fopt u(tindex, pricelevel, stocklevel) = uopt NEXT pricelevel NEXT stocklevel NEXT tindex REM Result tables PRINT "Optimal decision tables" PRINT "" PRINT "" FOR tindex = 0 TO tmax t = tindex * deltat PRINT " u(.) when time = "; t PRINT "" PRINT "Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 " PRINT "-------------------------------------------------------------" FOR stocklevel = 1 TO 15 PRINT USING "#####"; xlevel(stocklevel); PRINT " _"; FOR pricelevel = 1 TO 9 PRINT USING "######"; u(tindex, pricelevel, stocklevel); NEXT pricelevel PRINT "" NEXT stocklevel PRINT "" PRINT "" INPUT z NEXT tindex PRINT "" PRINT "Optimal result tables" PRINT "" PRINT "" FOR tindex = 0 TO tmax t = tindex * deltat PRINT " f(.) when time = "; t PRINT "" PRINT "Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 " PRINT "-------------------------------------------------------------" FOR stocklevel = 1 TO 15 PRINT USING "#####"; xlevel(stocklevel); PRINT " _"; FOR pricelevel = 1 TO 9 PRINT USING "######"; f(tindex, pricelevel, stocklevel); NEXT pricelevel PRINT "" NEXT stocklevel PRINT "" PRINT "" INPUT z NEXT tindex PRINT #2, "Optimal decision tables" PRINT #2, "" PRINT #2, "" FOR tindex = 0 TO tmax t = tindex * deltat PRINT #2, " u(.) when time = "; t PRINT #2, "" PRINT #2, "Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 " PRINT #2, "-------------------------------------------------------------" FOR stocklevel = 1 TO 15 PRINT #2, USING "#####"; xlevel(stocklevel); PRINT #2, " _"; FOR pricelevel = 1 TO 9 PRINT #2, USING "######"; u(tindex, pricelevel, stocklevel); NEXT pricelevel PRINT #2, "" NEXT stocklevel PRINT #2, "" PRINT #2, "" REM INPUT z NEXT tindex PRINT "" REM INPUT z PRINT "" PRINT #2, "Optimal result tables" PRINT #2, "" PRINT #2, "" FOR tindex = 0 TO tmax t = tindex * deltat PRINT #2, " f(.) when time = "; t PRINT #2, "" PRINT #2, "Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 " PRINT #2, "-------------------------------------------------------------" FOR stocklevel = 1 TO 15 PRINT #2, USING "#####"; xlevel(stocklevel); PRINT #2, " _"; FOR pricelevel = 1 TO 9 PRINT #2, USING "######"; f(tindex, pricelevel, stocklevel); NEXT pricelevel PRINT #2, "" NEXT stocklevel PRINT #2, "" PRINT #2, "" REM INPUT z NEXT tindex CLOSE #2 END ************************************************************************************ CCF harvest optimization via stochastic dynamic programming Peter Lohmander Example 2 2017-11-12 Optimal decision tables u(.) when time = 0 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 0 0 0 0 0 0 0 17 17 72 _ 0 0 0 0 0 0 25 42 42 107 _ 0 0 0 0 0 0 59 77 77 150 _ 0 0 0 0 0 43 103 120 120 199 _ 0 0 0 0 0 92 152 169 169 248 _ 0 0 0 0 0 141 201 218 218 291 _ 0 0 0 0 0 185 244 261 261 326 _ 0 0 0 0 0 220 279 296 296 352 _ 0 0 0 0 26 245 305 322 322 369 _ 0 0 0 0 43 263 322 339 339 381 _ 0 0 0 0 55 274 334 351 351 388 _ 0 0 0 0 62 282 341 358 358 393 _ 0 0 0 0 66 286 346 363 363 396 _ 0 0 0 0 69 289 348 366 366 397 _ 2 2 2 2 71 291 350 367 367 u(.) when time = 10 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 0 0 0 0 0 0 0 17 17 72 _ 0 0 0 0 0 0 25 42 42 107 _ 0 0 0 0 0 0 59 77 77 150 _ 0 0 0 0 0 43 103 120 120 199 _ 0 0 0 0 0 92 152 169 169 248 _ 0 0 0 0 0 141 201 218 218 291 _ 0 0 0 0 44 185 244 261 261 326 _ 0 0 0 0 78 220 279 296 296 352 _ 0 0 0 0 104 245 305 322 322 369 _ 0 0 0 0 121 263 322 339 339 381 _ 0 0 0 0 133 274 334 351 351 388 _ 0 0 0 0 140 282 341 358 358 393 _ 0 0 0 0 145 286 346 363 363 396 _ 0 0 0 0 148 289 348 366 366 397 _ 2 2 2 2 149 291 350 367 367 u(.) when time = 20 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 0 0 0 0 0 0 0 17 17 72 _ 0 0 0 0 0 0 25 42 42 107 _ 0 0 0 0 0 0 59 77 77 150 _ 0 0 0 0 0 43 103 120 120 199 _ 0 0 0 0 0 92 152 169 169 248 _ 0 0 0 0 0 141 201 218 218 291 _ 0 0 0 0 44 185 244 261 261 326 _ 0 0 0 0 78 220 279 296 296 352 _ 0 0 0 0 104 245 305 322 322 369 _ 0 0 0 0 121 263 322 339 339 381 _ 0 0 0 0 133 274 334 351 351 388 _ 0 0 0 0 140 282 341 358 358 393 _ 0 0 0 0 145 286 346 363 363 396 _ 0 0 0 0 148 289 348 366 366 397 _ 2 2 2 2 149 291 350 367 367 u(.) when time = 30 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 0 0 0 0 0 0 17 17 17 72 _ 0 0 0 0 0 0 42 42 42 107 _ 0 0 0 0 0 34 77 77 77 150 _ 0 0 0 0 0 78 120 120 120 199 _ 0 0 0 0 0 127 169 169 169 248 _ 0 0 0 0 49 176 218 218 218 291 _ 0 0 0 0 93 219 261 261 261 326 _ 0 0 0 0 128 254 296 296 296 352 _ 0 0 0 0 153 280 322 322 322 369 _ 0 0 0 0 171 297 339 339 339 381 _ 0 0 0 0 182 309 351 351 351 388 _ 0 0 0 0 189 316 358 358 358 393 _ 0 0 0 0 194 320 363 363 363 396 _ 0 0 0 0 197 323 366 366 366 397 _ 2 2 2 2 199 325 367 367 367 u(.) when time = 40 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 0 0 0 0 0 17 17 17 17 72 _ 0 0 0 0 0 42 42 42 42 107 _ 0 0 0 0 0 77 77 77 77 150 _ 0 0 0 0 43 120 120 120 120 199 _ 0 0 0 0 92 169 169 169 169 248 _ 0 0 0 0 141 218 218 218 218 291 _ 0 0 0 0 185 261 261 261 261 326 _ 0 0 0 0 220 296 296 296 296 352 _ 0 0 0 0 245 322 322 322 322 369 _ 0 0 0 0 263 339 339 339 339 381 _ 0 0 0 0 274 351 351 351 351 388 _ 0 0 0 7 282 358 358 358 358 393 _ 0 0 0 12 286 363 363 363 363 396 _ 0 0 0 15 289 366 366 366 366 397 _ 2 2 2 16 291 367 367 367 367 u(.) when time = 50 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 0 0 47 47 47 47 47 47 47 72 _ 0 0 72 72 72 72 72 72 72 107 _ 0 0 107 107 107 107 107 107 107 150 _ 0 0 150 150 150 150 150 150 150 199 _ 0 0 199 199 199 199 199 199 199 248 _ 0 0 248 248 248 248 248 248 248 291 _ 0 0 291 291 291 291 291 291 291 326 _ 0 0 326 326 326 326 326 326 326 352 _ 0 0 352 352 352 352 352 352 352 369 _ 0 0 369 369 369 369 369 369 369 381 _ 0 0 381 381 381 381 381 381 381 388 _ 0 0 388 388 388 388 388 388 388 393 _ 0 0 393 393 393 393 393 393 393 396 _ 0 0 396 396 396 396 396 396 396 397 _ 0 0 397 397 397 397 397 397 397 Optimal result tables f(.) when time = 0 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 1988 1988 1988 1988 1988 1988 1988 2075 2247 72 _ 2705 2705 2705 2705 2705 2705 2741 3078 3500 107 _ 3505 3505 3505 3505 3505 3505 3772 4453 5219 150 _ 4330 4330 4330 4330 4330 4370 5069 6183 7381 199 _ 5101 5101 5101 5101 5101 5349 6536 8139 9827 248 _ 5751 5751 5751 5751 5751 6330 8008 10101 12280 291 _ 6257 6257 6257 6257 6257 7202 9317 11847 14461 326 _ 6621 6621 6621 6621 6621 7899 10362 13240 16203 352 _ 6869 6869 6869 6869 6876 8409 11127 14260 17478 369 _ 7030 7030 7030 7030 7051 8759 11653 14961 18354 381 _ 7132 7132 7132 7132 7166 8989 11998 15421 18929 388 _ 7196 7196 7196 7196 7239 9136 12217 15714 19295 393 _ 7233 7233 7233 7233 7285 9228 12355 15897 19524 396 _ 7250 7250 7250 7250 7314 9284 12440 16010 19666 397 _ 7198 7215 7233 7250 7331 9319 12492 16080 19753 f(.) when time = 10 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 1426 1426 1426 1426 1426 1426 1426 1509 1636 72 _ 1993 1993 1993 1993 1993 1993 2039 2308 2621 107 _ 2637 2637 2637 2637 2637 2637 2880 3403 3970 150 _ 3306 3306 3306 3306 3306 3374 3937 4780 5668 199 _ 3933 3933 3933 3933 3933 4207 5133 6338 7589 248 _ 4471 4471 4471 4471 4471 5043 6332 7901 9515 291 _ 4890 4890 4890 4890 4891 5786 7398 9291 11228 326 _ 5192 5192 5192 5192 5226 6380 8250 10401 12596 352 _ 5397 5397 5397 5397 5472 6814 8874 11213 13597 369 _ 5531 5531 5531 5531 5641 7113 9302 11771 14285 381 _ 5616 5616 5616 5616 5751 7309 9583 12137 14736 388 _ 5669 5669 5669 5669 5822 7434 9762 12371 15024 393 _ 5700 5700 5700 5700 5866 7512 9874 12517 15204 396 _ 5715 5715 5715 5715 5893 7560 9943 12607 15315 397 _ 5681 5694 5706 5719 5910 7590 9986 12662 15383 f(.) when time = 20 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 988 988 988 988 988 988 988 1077 1171 72 _ 1421 1421 1421 1421 1421 1421 1483 1710 1941 107 _ 1929 1929 1929 1929 1929 1929 2162 2577 2998 150 _ 2463 2463 2463 2463 2463 2546 3016 3669 4327 199 _ 2972 2972 2972 2972 2972 3244 3983 4904 5830 248 _ 3411 3411 3411 3411 3411 3944 4952 6142 7338 291 _ 3754 3754 3754 3754 3794 4567 5814 7244 8679 326 _ 4000 4000 4000 4000 4100 5064 6502 8123 9749 352 _ 4168 4168 4168 4168 4324 5428 7006 8767 10533 369 _ 4278 4278 4278 4278 4478 5678 7352 9209 11071 381 _ 4347 4347 4347 4347 4579 5842 7579 9499 11425 388 _ 4391 4391 4391 4391 4643 5946 7724 9684 11650 393 _ 4417 4417 4417 4417 4683 6012 7815 9800 11791 396 _ 4430 4430 4430 4430 4708 6052 7871 9872 11878 397 _ 4407 4417 4426 4436 4723 6077 7905 9916 11931 f(.) when time = 30 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 655 655 655 655 655 655 688 758 827 72 _ 969 969 969 969 969 969 1085 1257 1429 107 _ 1347 1347 1347 1347 1347 1375 1630 1942 2253 150 _ 1759 1759 1759 1759 1759 1885 2316 2803 3291 199 _ 2161 2161 2161 2161 2161 2461 3092 3778 4464 248 _ 2510 2510 2510 2510 2540 3040 3869 4755 5641 291 _ 2784 2784 2784 2784 2877 3554 4561 5624 6687 326 _ 2981 2981 2981 2981 3146 3965 5114 6318 7523 352 _ 3116 3116 3116 3116 3343 4266 5518 6826 8135 369 _ 3204 3204 3204 3204 3479 4472 5796 7175 8555 381 _ 3259 3259 3259 3259 3568 4608 5978 7404 8831 388 _ 3294 3294 3294 3294 3624 4694 6094 7550 9006 393 _ 3316 3316 3316 3316 3660 4748 6167 7642 9116 396 _ 3327 3327 3327 3327 3681 4781 6212 7698 9184 397 _ 3312 3320 3327 3334 3695 4802 6239 7733 9226 f(.) when time = 40 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 411 411 411 411 411 432 484 535 587 72 _ 609 609 609 609 609 674 801 928 1055 107 _ 858 858 858 858 858 1005 1236 1466 1697 150 _ 1139 1139 1139 1139 1144 1422 1783 2144 2505 199 _ 1422 1422 1422 1422 1468 1893 2401 2910 3418 248 _ 1673 1673 1673 1673 1793 2366 3022 3678 4334 291 _ 1873 1873 1873 1873 2083 2786 3574 4361 5149 326 _ 2020 2020 2020 2020 2313 3122 4015 4907 5799 352 _ 2121 2121 2121 2121 2482 3368 4337 5307 6276 369 _ 2187 2187 2187 2187 2598 3537 4559 5581 6603 381 _ 2229 2229 2229 2229 2675 3648 4704 5761 6817 388 _ 2256 2256 2256 2256 2723 3718 4797 5876 6954 393 _ 2272 2272 2272 2272 2754 3762 4855 5947 7040 396 _ 2282 2282 2282 2283 2772 3790 4891 5992 7093 397 _ 2273 2278 2283 2289 2784 3807 4913 6019 7125 f(.) when time = 50 Stock p1 p2 p3 p4 p5 p6 p7 p8 p9 ------------------------------------------------------------- 47 _ 0 0 64 169 274 380 485 590 695 72 _ 0 0 92 253 414 575 736 898 1059 107 _ 0 0 130 368 606 844 1082 1320 1557 150 _ 0 0 178 513 847 1182 1516 1850 2185 199 _ 0 0 233 676 1120 1563 2007 2450 2894 248 _ 0 0 288 841 1394 1947 2500 3053 3605 291 _ 0 0 336 987 1637 2287 2938 3588 4238 326 _ 0 0 375 1103 1831 2559 3287 4015 4744 352 _ 0 0 404 1189 1974 2759 3544 4328 5113 369 _ 0 0 423 1247 2071 2895 3719 4543 5367 381 _ 0 0 436 1286 2135 2985 3835 4684 5534 388 _ 0 0 444 1310 2176 3042 3908 4774 5640 393 _ 0 0 449 1326 2202 3078 3954 4831 5707 396 _ 0 0 452 1335 2218 3100 3983 4865 5748 397 _ 0 0 454 1341 2227 3114 4000 4887 5773