ClearAll["Global`*"]

evalDate=DateObject[{2023,12,31}]; (*datum vyhotnocení T*)
{k,p}={0.5,0.5}; (*parametry rozdělení negativně binomického pro smlouvu s jednotkovou expozicí*)
(*Počátek platnosti smlouvy je mezi 1.1.2015 a 31.12.2022, smouva trvá 1, 2 nebo 3 roky. Funkce fce1 generuje datum počátku a konce platnosti pojistné smlouvy.*)
fce1[]:=With[{p=RandomInteger[{0,(DateObject[{2022,12,31}]-DateObject[{2015,1,1}])[[1]]}]},{DateObject[{2015,1,p}],DateObject[{2015,1,p+RandomChoice[Range[1,3]]*365}] }]
(*Funkce fce2 generuje datum nastání a nahášení škody. *)
fce2[{p_,k_}]:=With[{nastani=RandomChoice[Range[1,(Min[evalDate,k]-p)[[1]]]],zpozdeni=Ceiling@RandomVariate[ExponentialDistribution[1/(2*365)]]},{DayPlus[p,nastani],DayPlus[p,nastani+zpozdeni]}]
(*Funkce log1, log2, log3 slouží jako pomocné funkce v metodě maximální věrohodnosti.*)
log1[x_,y_]:=Log[PDF[ExponentialDistribution[s],x]/CDF[ExponentialDistribution[s],y]]
log2[x_,e_,q_]:=Log[PDF[NegativeBinomialDistribution[w,z/(z+e*q-z*e*q)],x]]
log3[x_,e_,q_]:=Log[PDF[PoissonDistribution[l*e*q],x]]

Simulace údajů

SeedRandom[2023];
pocetSmluv=1000;
smlouvyPlatnost=Table[fce1[],pocetSmluv]; (*údaje o 1000 smlouvách, ve formátu ...,{počátek, konec},...*)
eI=Map[((Min[#[[2]]-#[[1]],evalDate-#[[1]]])[[1]])/365&,smlouvyPlatnost]; (*riziková expozice*)
pocetSkod=Map[RandomVariate[NegativeBinomialDistribution[k,p/(p+#-p*#)]]&,eI]; (*celkový počet škod na smlouvě*)
skody=Map[Table[fce2[smlouvyPlatnost[[#]]],pocetSkod[[#]]]&,Range[1,Length@eI]]; (*údje o škodách ve formátu ...,{...,{nastání, nahlášení},...},... první úrovně seznamu odpovídá smlouvám, druhá úroveň samotným škodám, prazdný seznam odpovídá smlově s 0 škodami*)
pocetNahlasenych=Map[Count[#,{x_,y_}/;y<evalDate]&,skody];
(*vybrané ty škody, které byly nahlášeny před datem evaluace*)
nahlaseneSkody=Select[Flatten[skody,1],#[[2]]<evalDate&];

Aplikace modelu

reportLag=Map[(#[[2]]-#[[1]])[[1]]&,nahlaseneSkody]; (*zpozdění nahlášení*)
bodUseknuti=Map[(evalDate-#[[1]])[[1]]&,nahlaseneSkody]; (*rozdíl data vyhodncení a data nastání*)
(*logaritmická věrohodnostní funkce pro odhad parametru rozdělení doby zpoždění, předpoklad: rozdělení doby zpoždění je exponenciálním*)
logFce1=Total@Map[log1[reportLag[[#]],bodUseknuti[[#]]]&,Range[Length@reportLag]];
Plot[logFce1,{s,0,0.005}] (*graf, ze kterého určíme počáteční odhad pro numerickou metodu hledání maxima*)
(*Maximálně věrohodný odhad pro rozdělení doby zpoždění*)
sOdhad=FindArgMax[{logFce1,s>0  &&s∈Reals},{s,0.0015}] 
(*Funkce integrály určuje hodnotu dvojného integrálu, v textu práce rovnost (2.1)*)
integraly[V_,Z_]=Integrate[Integrate[PDF[ExponentialDistribution[sOdhad[[1]]],y]/V,{y,0,Z-w},Assumptions->Z-w∈Reals],{w,0,V},Assumptions->V∈Reals];
(*prravděpodobnost, že škoda byla nahlášena*)
qI=Map[integraly[(Min[#[[2]]-#[[1]],evalDate-#[[1]]])[[1]],(evalDate-#[[1]])[[1]]]&,smlouvyPlatnost];
(*logaritmická věrohodnostní funkce pro odhad parametrů rozdělení počtu nahlášených škod, předpoklad: rozdělení počtu nahlášených škod je negativně binomické*)
logFce2=Total@Map[log2[pocetNahlasenych[[#]],eI[[#]],qI[[#]]]&,Range[pocetSmluv]];
(*Maximálně věrohodné odhady pro rozdělení počtu nahlášených škod*)
{odhadK,odhadP}=ArgMax[{logFce2,0<w&& 0<z<1 &&  w∈Reals && z∈Reals},{w,z}]
(*logaritmická věrohodnostní funkce pro odhad parametru rozdělení počtu nahlášených škod, předpoklad: rozdělení počtu nahlášených škod je Poissonovo*)
logFce3=Total@Map[log3[pocetNahlasenych[[#]],eI[[#]],qI[[#]]]&,Range[pocetSmluv]];
(*Maximálně věrohodný odhad pro rozdělení počtu nahlášených škod*)
odhadLambda=ArgMax[{logFce3,l>0 && l∈Reals},l]
(*V textu práce funkce (OverBar[Q])(j), počet náhlášených škod se řídí negativně binomickým rozdělením.*)
prumerNB[w_]:=(1/pocetSmluv)Total[Map[SetPrecision[PDF[NegativeBinomialDistribution[odhadK,odhadP/(odhadP+eI[[#]]*qI[[#]]-odhadP*eI[[#]]*qI[[#]])],w],3]&,Range[pocetSmluv]]]
(* (OverBar[Q])(j), počet náhlášených škod se řídí Poissonovým rozdělením.*)
prumerPo[w_]:=(1/pocetSmluv)Total[Map[SetPrecision[PDF[PoissonDistribution[odhadLambda*eI[[#]]*qI[[#]]],w],3]&,Range[pocetSmluv]]]
(*Relativní četnost počtu škod w*)
A[w_]:=SetPrecision[Count[pocetSkod,w]/pocetSmluv,3]
Module[{a=Range[0,10]},MapThread[Prepend,{{a,Map[A,a],Thread[prumerNB@a],Thread[prumerPo@a]},{"počet","pozorované","NBi","Po"}}]]//MatrixForm

(*V řádcích jsou hodnoty T(j)/(1-T(j)), (OverBar[Q])(j)/(1- (OverBar[Q])(j)) postupně pro Poissonův a negativně binomický případ.*)
Module[{a=Range@Max[pocetSkod], b=prumerNB[0],c=prumerPo[0],d=A[0]},MapThread[Prepend,{{a,Map[A,a]/(1-d),Thread[prumerNB@a]/(1-b),Thread[prumerPo@a]/(1-c)},{"počet","pozorované","NBi","Po"}}]]//MatrixForm

(*dodělat rozdělení nenahlášených*)

P=N[Mean[pocetSkod/eI]]

(odhadK(1-odhadP))/odhadP