Definice parametrů, systému a neznámých pro soustavu rovnic

(* škodní frekvence, parametr Poissonova rozdělení *)
lambda = 0.1;
(* parametry gamma rozdělení výše škod *)
xi = 10; nu = 2;
(* počet škod již nahlášených v tomto roce *)
m=0;
(* roční úroková míra *)
i=0.1;
(* pojistné v základní bonus-malus třídě *)
baseP = 1;
(* počet tříd, o které se pojistník při nahlášení nehody posune výše *)
a=3; 
(* relativní pojistná v systému *)
SYSTEM={0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2};
(* diskontní faktor *)
d = 1/(1+i);
(* neznámé pro výpočet současné hodnoty *)
pvVars = Table[Symbol["$pv"<>ToString@(i)],{i,Length[SYSTEM]}];
(* neznámé pro výpočet retenční hranice *)
rlVars = Table[Symbol["$rl"<>ToString@(i)],{i,Length[SYSTEM]}];

Definice funkcí potřebných pro výpočet

(* přechodová funkce *)
levelFunction[level_, k_] := If[k==0, Max[level - 1, 1], Min[level + a*k, Length[SYSTEM]]]

(* hustota gamma rozdělení výše škod *)
gammaDensity[y_] := 1/Gamma[nu]*(nu/xi)^nu * y^(nu - 1)*Exp[-y*nu/xi];

(* pravděpodobnost nenahlášení škody při dané retenční hranici *)
probNot[rl_] := Integrate[gammaDensity[y], {y, 0, rl}]

(* pravděpodobnost nahlášení právě k škod pro danou úroveň v systému při dané retenční hranici *)
probK[level_, k_, rl_] :=Sum[Exp[-lambda]*lambda^h/Factorial[h]*Binomial[h, k]*((1-probNot[rl])^k)*(probNot[rl]^(h-k)),{h, k,Infinity}];

(* průměrná frekvence nahlášených škod *)
avgLambda[level_, rl_] := Sum[k*probK[level, k, rl], {k, 0, Infinity}];

(* očekávaná cena nenahlášené nehody pro danou retenční hranici *)
ECost[rl_] := If[probNot[rl] > 0,(1/probNot[rl])*Integrate[y*gammaDensity[y], {y, 0, rl}], 0];

(* výpočet současné hodnoty pro danou úroveň a retenční hranici *)
presentValue[level_, rl_] := Module[{return},
If[rl ≠ 0,
return = SYSTEM[[level]]*baseP + 
d^(1/2)*ECost[rl]*(lambda-avgLambda[level, rl])+
d*Sum[probK[level, k, rl]*pvVars[[levelFunction[level, k]]], {k, 0, Quotient[Length[SYSTEM], a]}];
return += d*(1-Sum[probK[level, k, rl], {k, 0, Quotient[Length[SYSTEM], a]}])*pvVars[[Quotient[Length[SYSTEM], a] + 1]],
return = SYSTEM[[level]]*baseP + 
d*Sum[Exp[-lambda]*lambda^k/Factorial[k]*pvVars[[levelFunction[level, k]]], {k, 0, Quotient[Length[SYSTEM], a]}];
return += d*(1-Sum[Exp[-lambda]*lambda^k/Factorial[k], {k, 0, Quotient[Length[SYSTEM], a]}])*pvVars[[Quotient[Length[SYSTEM], a] + 1]]
];
return
];

(* výpočet retenční hranice pro danou úroveň a současné hodnoty *)
retention[level_, presentValues_] := Module[{levF1, levF0}, 
levF1 = levelFunction[level,m+1];
levF0 = levelFunction[level,m];
presentValues[[levF1]] - presentValues[[levF0]]
]

(* výpočet optimální retence pro systém s danými parametry *)
calculation[lambda_, xi_, nu_, m_, d_, a_, baseP_, SYSTEM_] := Module[{rls, rlsBefore, pvsSolved, pvs, i},
i = 0;
rls = Table[0, {i, Length[SYSTEM]}];
rlsBefore = Table[Null, {i, Length[SYSTEM]}];
While[
Not[SameQ[Round[rls, 0.00001], Round[rlsBefore, 0.00001]]],
rlsBefore = rls;
pvsSolved = Solve[MapThread[Equal,
{Table[presentValue[level, rls[[level]]], {level, Length[SYSTEM]}],pvVars}],
pvVars
];
pvs = pvVars /.pvsSolved[[1]];
rls = Table[retention[level, pvs], {level, Length[SYSTEM]}];
Print["Iteration ", i, ": ", rls];
i += 1;
];
Print["Completed, optimal retentions are: ", rls]
]

Výpočet optimální retence

calculation[lambda, xi, nu, m, d, a, baseP, SYSTEM]