EconBase
← Back to paper

Infinitely Stochastic Micro Forecasting

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

85,231 characters · 19 sections · 76 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Infinitely Stochastic Micro Forecasting

abstractForecasting costs is now a front burner in empirical economics. We propose an unconventional tool for stochastic prediction of future expenses based on the individual (micro) developments of recorded events. Consider a firm, enterprise, institution, or state, which possesses knowledge about particular historical events. For each event, there is a series of several related subevents: payments or losses spread over time, which all leads to an infinitely stochastic process at the end. Nevertheless, the issue is that some already occurred events do not have to be necessarily reported. The aim lies in forecasting future subevent flows coming from already reported, occurred but not reported, and yet not occurred events. Our methodology is illustrated on quantitative risk assessment, however, it can be applied to other areas such as startups, epidemics, war damages, advertising and commercials, digital payments, or drug prescription as manifested in the paper. As a theoretical contribution, inference for infinitely stochastic processes is developed. In particular, a non-homogeneous Poisson process with non-homogeneous Poisson processes as marks is used, which includes for instance the Cox process as a special case.
smallKeywords: stochastic prediction, infinitely stochastic process, marked process, time-varying models, dynamic panel data, resampling, risk valuation

Introduction

Human as well as monetary losses and uncertainty about their extent are one of the main sources of risk. A probabilistic prediction of the future monetary losses lies on the cutting edge of quantitative risk assessment, for instance, valuation of operational risk in banking or reserving risk in insurance, while the study of human losses is of particular interest in conflict solution and epidemics modeling. We propose a general prediction methodology, which together with the underlying stochastic procedures is applicable to various areas as demonstrated by case examples later on.

Let us define the general structure of our model on the basis of a financial example. The event's `lifetime' can be described as follows: The $i$th loss occurs at the occurrence time, which is denoted by $T_i$. Such a loss is often reported (e.g., to a financial company) not immediately after the event, but for various reasons, after the reporting delay (waiting time) $W_i$, which is the time difference between the occurrence epoch (event time) and the observation epoch (reporting time). Furthermore, $Z_i=T_i+W_i$ stands for the $i$th reporting (notification) time. The contemplated cash flows are visualized in Figure (ref), which elucidates the whole framework behind the loss reporting process together with the time developments of the losses.

figure[figure omitted — 3,995 chars of source]

Our observation history for the reported losses is a time interval $[0,a]$, where $a$ is the present time. The main aim is to predict the losses, which are going to be reported in the future time horizon $(a,b]$, and simultaneously to predict the development of the losses within the time interval $(a,b]$, that have already occurred before time $a$, but are not settled yet. Some of them are already incurred (i.e., occurred within $[0,a]$, but will be reported after time point $a$). The observed loss data are truncated in the way that we observe only the reported losses, i.e., $Z_i\leq a$. Without loss of generality, we assume that the reporting times are chronologically ordered such that $Z_{i_1}\leq Z_{i_2}$ for $i_1<i_2$. After some internal process is carried out, the company pays $N_i(t)$ payments till time point $t$, or $N_i(\infty)$ payments in order to fully settle the $i$th loss. The amount of the $k$th payment within the $i$th loss paid at time $U_{i,k}$ is represented by $X_{i,k}$, for $k=1,\ldots,N_i(a)$. The time window from the reporting time $Z_i$ up to the last observed (available) time $a$ for the $i$th loss has length $V_i$, i.e., $V_i=a-Z_i$. One can think of the reporting times $Z_i$'s as the arrival times of the counting process $\{M(t)\}_{t\geq 0}$ and the payment times $U_{i,k}$'s as the arrival times of the counting processes $\{N_i(t)\}_{t\geq 0}$ for $i\in\mathbb{N}$. Assuming that there are $i=1,\ldots,M(a)$ losses already reported, we observe a collection \[ \{T_i,Z_i,\{U_{i,j}\}_{j=1,\ldots,N_i(a)},\{X_{i,j}\}_{j=1,\ldots,N_i(a)}\}_{i=1,\ldots,M(a)} \] or, alternatively and equivalently, $\{Z_i,W_i,\{N_i(t)\}_{t\in[0,a]},\{X_{i,j}\}_{j=1,\ldots,N_i(a)}\}_{i=1,\ldots,M(a)}$.

Motivation and applications

The proposed class of models---infinitely stochastic processes---is a very rich and general class that nests, for examples, doubly stochastic (Cox) processes. Our approach and results are motivated in the context of several applications taken from the empirical economics literature.

\paragraph{Case 1: Operational risk} Banks and other financial institutions have to face operational risk covering fraud, system failures, security, privacy protection, terrorism, legal risks, employee compensation claims, physical (e.g., infrastructure shutdown) or environmental risks. As pointed out in OR2007, some large banks prefer to use their own formal definition of operational risk. For instance, Deutsche Bank (2017) defines operational risk as “the risk of loss resulting from inadequate or failed internal processes, people and systems or from external events, and includes legal risk.” Recent developments for operational risk prediction---comprehensively summarized by BLM2018---reveal that challenges like truncated data and non-homogeneous processes have to be handled in operational risk modeling. The empirical literature, for instance Cohen2018, suggests that operational risk capital models can be based on the loss distribution approach. Generally, a loss $i$ corresponding to operational risk is occurred at $T_i$, but is internally reported later at $Z_i$. Consequently, compensations need to carried out by the bank to the affected side. For the loss $i$, the compensations $X_{i,k}$'s are going to be paid at the times $U_{i,k}$'s. The bank is then required (e.g., by the Third Basel Accord) to quantify the future distribution of losses (measured through their compensations) belonging to operational risk.

\paragraph{Case 2: War damages} Modeling of the evolution of national and international conflicts is a long standing strand of research GleditschMetternichRuggeri2014, although data collection is a very challenging task, see Arnold2019 and Cressey2008. Based on the comprehensive databases as COPDAB Azar1980, MIDLOC Braithwaite2010, or PRIO Hallberg2012, main approaches still remain to be classic econometric linear ones with a list of exogenous factors or those based on the hazard models, cf. CollierHoefflerSoderbom2004, Schrodt2014, Clauset2014 and BakkeGreenhillWard2010, HarrisonWolf2012. As the list of current approaches was criticized by Schrodt2014, namely “garbage can models that ignore the effect of collinearity”, “complex models without understanding the underlying assumptions” or “linear statistical monoculture”, we believe that our model can make a useful step forward in the conflict prediction. Using the proposed model, one considers each point $T_i$ in the $M$ process as the beginning of the tension between regions/countries. The followed up $Z_i$ is the official beginning of the conflict, through the official notice or the first armed intrusion. This point (as the mark) starts a spread the armed conflict over a series of battles $k$ at time points $U_{i,k}$ that take $X_{i,k}$ lifes. The main assumptions of the process are fully in-line with the nature of the war: (ref) “there will always be another war” and (ref) “each war has an end”.

\paragraph{Case 3: Epidemics} Proper modeling and prediction of the spread of epidemic is a very important strand of literature, in particular in the view of recent H1N1 and Ebola epidemics. Classical models arising from modeling the online diffusions are those based on Susceptible-Infected-Recovered model introduced by Kermack27, LindaAllen2008, often extended to the stochastic case as in bobashev2007, PingYan2008, and others. Newly, these models were linked to the Hawkes processes in which spread of the disease in one population has been investigated, see SIRHawkes2018. Considering several populations (neighbor regions, countries, flight connections, etc.), the proposed model is a natural flexible extension, in which each point $T_i$ of the process is the infection of an individual in the region $i$. The delay $W_i$ is thus the incubation period, after which the process of infection individuals in the population $i$ starts, with individuals being infected at the points $U_{i,k}$. The values $X_{i,k}$ may be considered as the severity of the illness also converted to monetary quantities.

\paragraph{Case 4: Drug prescription} Health care expenditures have become one of the most serious issues of the modern society and prescribed medicaments seem to form one of the fastest growing component of the health insurance expenses. Managed care organizations encourage physicians to be more cost-conscious, see, e.g., miller. They use financial incentives to induce physicians to reduce expenses while maintaining the quality of medical care. General practitioners (GP) comprise a significant part of the health care system and influence importantly the total amount of insurance money spent. Therefore, the prescription behaviour of GPs is of utmost importance. It has been studied from the point of view of the pharmaceutical firms in several studies, see gonul, manchanda, and other references therein; or from the point of view of the health insurance companies, e.g., HPH2017. Furthermore, the prescription patterns and the influencing factors have been analyzed in ekedahl, rokstad, watkins, or Caldbick2015. One may think of a spread of disease/illness as the occurrence time $T_i$ of event and, correspondingly, a visit to the GP as the reporting time $Z_i$. Expenses for the prescribed drugs---sometimes more than one medical examination by the GP is needed (at times $U_{i,k}$)---are then the event payments $X_{i,k}$. After all, the responsible health care financing organization is interested to know the future expenditures for prescribed medicaments by the GPs within a predetermined time horizon.

\paragraph{Case 5: Startups} A recent entrepreneurship bloom prompts for another straightforward application of the proposed micro forecasting method. Many well-known multinational companies leading the global market these days have begun their business in terms of small and locally based startups with only very limited human, social, and financial capital. These are, although, crucial factors for the future startup performance and its ability to survive BPTW2004. Especially the last one---the financial capital---turns out to play the most significant role for establishing an entrepreneur on the global market CG2008. The initial financial capital is, however, usually not sufficient to start operations at the desired scale as the credit constraints for bank loans are too strict HJR1994. Therefore, the startup team looks for additional external sources of equity capital (such as external support, collaboration, fund raising, etc.), which is credited later over time. When modeling the overall impact of the startup in terms of its frontier production, this additional capital should be also considered AKS1977. Thus, the future cash flow is of the main interest. Using our terminology, a new entrepreneur $i$ starts with its business after some waiting time $W_i$ and additional capital amounts $X_{i,k}$'s arrive at the times $U_{i,k}$'s. The whole scenario can be analogously also adapted, for instance, for modeling the human capital, social investments, or business performance of the startups.

\paragraph{Case 6: Advertising and commercials} In recent years, television and internet advertising have become increasingly tailored to individuals. Television commercials simply rely on so-called contextual advertising, where ads are chosen based upon the broadcast contents LTW2015. More sophisticated ad placement techniques use online behavioral advertising or targeting advertising, which are typical for internet commercials and are based on the browsing history, online activities, or web searches GT2011. Shopping behavior and shopping patterns have recently been analyzed by economists in order to better understand the process through which consumers search for their preferred options BURDA2012,XIAO2018. Let us concentrate on the perspective of a company selling a specific product over the Internet. A moment, when some ad is placed over the Internet by another company running some website, can be thought of the occurrence time $T_i$. The owner of the website is directly or indirectly paid for the advertisement by the product-selling company. Then, the reporting time $Z_i$ naturally corresponds to the time when the ad is displayed or when it is recognized by some user (for example, the first click on it). If the user proceeds with a purchase in the online store advertised by the ad, the payment time and the payment amount represent the mark of an event development process. The product-selling company can consequently predict their future income based on the advertising and, thus, judge the efficiency of their commercial product placement.

\paragraph{Case 7: Apple Pay} Digital payment platforms are multi-sided and layered modular artifacts that primarily mediate payment transactions between payers and payees Kazan2015. Apple Pay as a payment method has become increasingly popular in everyday life during recent years. As remarked by LIU2019268, it simply boosts satisfaction through elevated coolness in a successful encounter. LIU2015372 examined recent changes in the payment sector in financial services, specifically related to mobile payments that enable new channels for consumer payments for goods and services purchases, and other forms of economic exchange. Although, Apple Pay serves merely as a proxy and mediator between cardholders' (card issuer) and merchants' (acquirer) bank accounts, a card issuer (e.g., a bank) is obliged to pay a portion from the payment for such a service. Therefore, from the card issuer's perspective, it is of interest to determine the amount of charges for the Apple Pay service during the future time period. Here, the date of issuing a payment card can be considered as the occurrence time $T_i$ of event, the date of registering the Apple Pay is the reporting time $Z_i$, and the corresponding purchases are the payment amounts $X_{i,k}$ realized at payment times $U_{i,k}$.

\paragraph{Case 8: Actuarial claims reserving} An insurance company needs to predict future claims with corresponding payments and, additionally, future payments coming from already occurred claims, which do not have to be necessarily reported, and this is due to the current regulatory framework for insurance supervision (e.g., Solvency II Directive or Swiss Solvency Test). A lifetime of a claim can be characterized by the following variables that are driving the claim process: A claim $i$ (i.e., a loss) occurs at the accident time $T_i$, however the insurance company is notified with some delay (heavy injuries that did not allow insured person report the claim, too light damages of the vehicle that allow an insured person to postpone a report, etc.) at the reporting time $Z_i$; the corresponding claim payments $X_{i,k}$'s are going to be paid by the insurance company to the insured at the payment times $U_{i,k}$'s. Finally, the distribution of cumulative payments within a predetermined time window has to be predicted in order to settle the required claim reserves (e.g., Value at Risk or Expected Shortfall at $99.5\%$). Below, we concentrate in more details on the claims reserving task and exemplify the proposed methodology through analyses of two insurance lines of business in order to demonstrate practical efficiency of our prediction method.

Outline

This paper is structured as follow: Next section introduces the data and main stochastic objects we intend to model. Section (ref) provides the assumptions and theory for occurrence and reporting times; number of payments; reporting delays; and payment amounts in different subsections. Section (ref) contains a practical application to the actuarial data. Afterwards, our conclusion follows. The proofs of our theoretical results are put in the Appendix.

Granular loss reserving

A classical actuarial problem called claims reserving is elaborated from an emerging perspective. In contrast to the traditional claims reserving techniques based on aggregated information from historical data, our approach relies on granular individual claim-by-claim data and contributes to increase in the prediction's precision.

Claims (loss) reserving in insurance determines a sufficient amount of money, that needs to be put aside from the premium, to cover future claim (loss) payments. The main issue is to estimate/predict these claims reserves, which should be held by the insurer in order to meet all future claims arising from policies currently in force and policies written in the past. Claims reserving is a classical problem in non-life insurance, sometimes also called general insurance (in UK) or property and casual insurance (in USA). A non-life insurance policy is a contract between the insurer and the insured. The insurer receives a deterministic amount of money, known as premium, from the insured in order to obtain a financial coverage against well-specified randomly occurred events. If such an event (claim) occurs, the insurer is obliged to pay in respect of the claim a claim amount, also known as a loss amount. In layman's terms, if an accident happens to an insured person, he or she goes to the insurance company to request a claim payment. The insurance company pays this claim amount from the loss reserves. In many cases several payments are performed for a single accident, for instance, further health problems, hidden damages of the car that were not visible by the first inspection, etc.

Claims reserving methods based on aggregated data from so-called run-off triangles are predominantly used to calculate the claims reserves, see EV2002 or wutrich_kniha for an overview. Such models are not based on the particular claims or accidents, but rather on the aggregated overall payments through some predefined period, typically one year. These conventional reserving techniques have series of disadvantages: loss of information from the policy and the claim's development due to the aggregation; usually small number of observations in the aggregated data; only few observations for recent accident years; various assumptions of independence, which can sometimes be unrealistic or at least questionable; and sensitivity to the most recent paid claims, see PH2013 or PO2014 for some recent developments. In order to overcome the above mentioned deficiencies or imperfections, micro (granular) loss reserving methods for individual claim-by-claim data need to be derived. Moreover, estimation of the whole distribution of the total future payments is a crucial part of the risk valuation process.

Current status

To estimate the distribution of reserves means to predict future cash flows and their uncertainty. On the top of that, this is becoming compulsory by the law due to the introduction of new supervisory guidelines.

The loss reserving approaches based on individual/micro-level/granular/claim-by-claim data do not represent the mainstream in the reserving field. First attempts within the reserving framework of incorporating the claim information for reporting delays were using a Bayesian approach Jewell1989,Jewell1990 or an empirical-Bayes approach WTC1984. Substantial branch of the individual loss reserving methods, that are based on a position dependent marked Poisson process, involves work of Arjas, Norberg,Norberg1999, and HaastrupArjas. A Markov model for granular loss reserving was proposed by Hesselager1994. Including claims features to specify the model components within the setup of the marked point processes was revisited by Larsen. Empirical investigation by AntonioPlat indicates that individual reserving provides a better accuracy compared to some selected aggregated models. A discrete time formulation instead of the continuous time point process description was suggested by GodecharleAntonio and PigeonAntonioDenuit. Besides that, ZhaoZhouWang and ZhaoZhou proposed semiparametric techniques from survival analysis. Several case studies of the individual reserving approaches can be found in TaylorMcGuireSullivan. Machine learning techniques in the individual claims reserving were elaborated by W2016. Furthermore, VW2016 pointed out a gain in using the individual methods by employing non-stationarity. Cox processes were utilized by BLT2016.

Practical loss prediction techniques often forfeit diagnostics of the theoretical models' assumptions. This is overcome by our approach, where we deal with nonlinear continuous time Markov environments HS2009. Generally, times of events together with accompanying measures can be analyzed as marked point processes, e.g., in case of ultra-high-frequency data, see Engle2000. Stochastic methods for modeling the total claim amount via marked Poisson cluster models have been recently proposed by BWZ2018. Here, the marks can take values only in a finite-dimensional space and have a common distribution. Our approach allows for point processes as marks, which is very suitable for a practical application of unrestricted number of payments, and enables time-varying distributions for, e.g., payment dates, reporting delays, or payment amounts.

Main goals

This paper contributes to the literature by aiming at using all the available information in the data. The proposed model motivated by the data controls for dependencies (between different payment amounts, between payments amounts and reporting delays, between reporting delay and accident date, etc.) in a simple and natural way. We assume a marked non-stationary Poisson process for the time ordered reporting dates; flexibly parametrized conditional distribution of the reporting delay and payment sizes given the accident date; or a non-homogeneous Poisson process in the role of a process' mark for the number of payments. All the models in this paper are supported with the asymptotic theory, and are combined in an omnibus model. Application to the true data strongly outperforms classical models in both point and interval forecast of the reserved losses. Up to the best of our knowledge, this is the first time where all the possible cross and temporal dependencies of the claim data are taken into account.

Data

Nowadays, modern databases and computer facilities provide a foundation for loss reserving based on individual data. There is no more reason to rely on the reserving techniques using aggregated data only. We possess the unique database from the Guarantee Fund of the Czech Insurer's Bureau for car insurance which consist of claims developments from the beginning of 2004 up to the end of 2016. Each record in the data set contains:

itemize• Claim ID (if one claim is associated with more payments, each payment has a separate entry); • Type of claim, which can be either bodily injury or material damage; • Accident time (occurrence); • Reporting time (notification); • Date of payment, when the payment is credited to the client's bank account; • Amount of payment.

All in all we have 4450 claims comprising of 10820 payments for bodily injuries and 30545 claims distributed into 35642 payments for material injuries within the investigated time interval. For back-testing purposes, we only use the data up to the end of 2015 to construct the prediction. The data from 2016 are only employed for comparison purposes with the obtained results.

Theoretical framework

The size of the claims reserves protects insurance company against the future losses. Bearing this practical issue in mind, we propose a series of theoretical models describing each component of the claim's chain. We assume that the reporting dates $Z_i$'s follow a non-homogeneous Poisson process with a parametric intensity function; the reporting delays $W_i$'s follows a time-varying continuous parametric distribution conditional on the reporting dates; the payment dates for each claim $i$ are represented by arrival times of a non-homogeneous Poisson process $N_i(t)$, where $N_i$ is a mark of $M$; and the payment amounts $X_{i,j}$'s are modeled similarly to the reporting delays via a time-varying parametric conditional distribution. All the models are then brought together under one umbrella in the empirical study.

Recently, GIESECKE2018 have discussed marked point processes with applications in finance and economics to model the timing of defaults, corporate bankruptcies, market transactions, unemployment spells, births, and mortgage delinquencies. They developed likelihood estimators for the parameters of a marked point process and incompletely observed explanatory factors that influence the arrival intensity and mark distribution, although they presumed only finite dimensional marks. We go beyond therein investigated models by assuming infinite dimensional marks via different stochastic framework. We establish an approximation to the likelihood and analyze the convergence and large-sample properties of the associated estimators. Numerical results illustrate the behavior of our estimators.

Recall that our primary practical goal is to model and, consequently, to simulate a distribution of the sum of future payments within the time period $(a,b]$. Besides that, the secondary practical goal is to back-predict claims that have already occurred, but are still not reported. As a theoretical by-product, we investigate marked non-homogeneous Poisson processes with infinite dimensional marks, which has not been done yet.

Occurrence and reporting times

As for the insurance company it is not so relevant, when the accident has happened, but rather when it has been reported, as on that date the whole procedure of the claim payments starts.

We proceed to the assumptions on the reporting dates $Z_i$'s that are needed for showing existence, uniqueness, and asymptotic properties of the proposed estimators. Taking into account, that different claims are supposed to be unrelated, it is natural to assume that the time differences $\{Z_i-Z_{i-1}\}_{i\in\mathbb{N}}$ are independent ($Z_0\equiv 0$). Thus, the reporting dates $Z_i$'s can be viewed as the arrival times of a counting process with independent increments. A reasonable and parsimonious representative would be a non-homogeneous Poisson process.

assumpMThe time ordered reporting times $\{Z_i\}_{i\in\mathbb{N}}$ are arrival times of a non-homoge\-neous Poisson process $\{M(t)\}_{t\geq 0}$ with a parametric intensity $\psi(t;\bm{\rho})>0$ such that $M(t)=\sum_{i=1}^{\infty}\mathbbm{1}\{Z_i\leq t\}$, $\bm{\rho}\in{\boldsymbol R}\subseteq\mathbb{R}^q$, and ${\boldsymbol R}$ is an open convex set.

The reporting epochs $\{Z_i\}_{i\in\mathbb{N}}$ are reversely determined by the counts $\{M(t)\}_{t\geq 0}$ such that $Z_i=\inf_{t\geq 0}\{M(t)\geq i\}$. The intensity $\psi(t;\bm{\rho})$ can be considered as a risk exposure for an accident reporting (not occurring) in time $t$. Although, one may still argue that it should be more convenient to assume that the occurrence (accident) times, and not the reporting times, should form the arrival times of some non-homogeneous Poisson process. This is indeed in concordance with the parametric time-varying conditional density $f_W$ for the reporting delay $W_i$ (given $Z_i=z$) defined later on, which results in the fact that $\{T_i\}_{i\in\mathbb{N}}$ are arrival times of another non-homogeneous Poisson process having the intensity

equation[equation omitted — 146 chars of source]

because of the displacement theorem Kingman1993. Thus, $\mu(t;\bm{\rho},\bm{\vartheta})$ is just a risk exposure for an accident occurring in time $t$.

We should emphasize that similarly to the below derived statistical inference for the non-homogeneous Poisson process, many other authors dealt with consistent estimation of the process intensity. To mention at least some of them, we refer to Konecny1987, SCHOENBERG2005, Waag2007, WG2009, coeurjolly2014, and PDV2017, although sometimes in a more general setup. There are two main reasons why we derive consistency and asymptotic normality of the intensity estimator in a different fashion: First, it is of a practical interest to require simple assumptions, which are easily verifiable and allowing for a huge class of parametric intensities. Second, our theoretical results and ways of proving them serve as an intermediate product for developing a suitable statistical inference for the marked non-homogeneous Poisson process with marks being non-homogeneous Poisson processes (discussed in Subsection (ref)).

Since we consider a fully parametric approach, it is firstly necessary to estimate the unknown parameter $\bm{\rho}$. We employ the maximum likelihood (ML) approach for the arrival times. The unconditional likelihood in case of $M(t)$, when the last observable (deterministic) time is $t$, has the form

equation[equation omitted — 167 chars of source]

where $\Psi(t;\bm{\rho}):=\int_{0}^t\psi(z;\bm{\rho})\mbox{d} z=\mathsf{E} M(t)$ is a cumulative intensity function.

Maximizing the log-likelihood function

equation[equation omitted — 132 chars of source]

with respect to $\bm{\rho}$ provides an ML estimator $\widehat{\bm{\rho}}$. The true value $\bm{\rho}_0\in{\boldsymbol R}$ of the unknown parameter $\bm{\rho}$ is supposed to uniquely maximize $\mathsf{E}_{\bm{\rho}}\ell(\bm{\rho};\bm{Z},t)$. The uniqueness of $\bm{\rho}_0$ is essential for identifiability and, consequently, for consistency and asymptotic normality of the ML estimator White1982. These properties are proved later on. Although, one has to realize that we are dealing with not independent and not identically distributed (n.i.n.i.d.) random variables (i.e., arrival times). Let us define \[ h(Z_i;\bm{\rho},t):=\frac{1}{M(t)}\Psi(t;\bm{\rho})-\log \psi(Z_i;\bm{\rho}), \] which means that

equation[equation omitted — 129 chars of source]

Assume the following to hold with respect to the $h(Z_i;\bm{\rho},t)$ functions in order to obtain a sensible estimator (i.e., consistent and asymptotically normal).

assumpM$h(z;\bm{\rho},t)$ is convex in $\bm{\rho}\in{\boldsymbol R}$ for all $0<z<t$.

The convexity of $h$ from Assumption (ref) contributes to assurance that there exists a unique optimum, namely some function $\widehat{\bm{\rho}}\equiv\widehat{\bm{\rho}}(t)$.

For simplicity of further notations, let us denote $[\cdot][\cdot]^{\top}\equiv[\cdot]^{\otimes 2}$, $\partial_{\bm{\rho}_0}\equiv\frac{\partial}{\partial\bm{\rho}}\left[\cdot\right]_{\bm{\rho}=\bm{\rho}_0}$, $\partial^2_{\bm{\rho}_0}\equiv\frac{\partial^2}{\partial\bm{\rho}\partial\bm{\rho}^{\top}}\left[\cdot\right]_{\bm{\rho}=\bm{\rho}_0}$, $\partial_{\bm{\rho}_0,i}\equiv\frac{\partial}{\partial\rho_i}\left[\cdot\right]_{\bm{\rho}=\bm{\rho}_0}$, and $\partial^2_{\bm{\rho}_0,i,j}\equiv\frac{\partial^2}{\partial\rho_i\partial\rho_j}\left[\cdot\right]_{\bm{\rho}=\bm{\rho}_0}$ for $\bm{\rho}=[\rho_1,\ldots,\rho_q]^{\top}$. Symbol $\bm{0}$ stands for a zero vector and ${\boldsymbol I}$ means an identity matrix with a suitable dimension. Furthermore for $t>0$, let us define an information matrix $\mathcal{I}(t;\bm{\rho}_0):=\int_0^t\frac{\{\partial_{\bm{\rho}_0}\psi(z;\bm{\rho})\}^{\otimes 2}}{\psi(z;\bm{\rho}_0)}\mbox{d} z$ and a matrix $\mathcal{K}(t,\bm{\rho}_0):=\frac{1}{\sqrt{\psi(t;\bm{\rho}_0)}}\left[\partial^2_{\bm{\rho}_0}\psi(t;\bm{\rho})-\frac{\{\partial_{\bm{\rho}_0}\psi(t;\bm{\rho})\}^{\otimes 2}}{\psi(t;\bm{\rho}_0)}\right]$.

assumpM$\frac{\partial^2}{\partial\bm{\rho}\partial\bm{\rho}^{\top}}\psi(\cdot;\bm{\rho}):\,\mathbb{R}^+\to\mathbb{R}^{q\times q}$ is continuous for all $\bm{\rho}\in{\boldsymbol R}$ and there exist Lebesgue-integrable functions $m_{1,i}$ and $m_{2,i,j}$ such that \[ \left|\frac{\partial}{\partial\rho_i}\psi(t;\bm{\rho})\right|\leq m_{1,i}(t)\quad\mbox{and}\quad\left|\frac{\partial^2}{\partial\rho_i\partial\rho_j}\psi(t;\bm{\rho})\right|\leq m_{2,i,j}(t) \] for all $\bm{\rho}\in{\boldsymbol R}$, almost every $t>0$, and $i,j=1,\ldots,q$.

The definition of $\mathcal{I}(t;\bm{\rho}_0)$, the uniqueness of the true value $\bm{\rho}_0\in{\boldsymbol R}$, and the differentiability of $\psi(t;\cdot)$ from Assumption (ref) ensure that $\mathcal{I}(t;\bm{\rho}_0)$ is positive definite. Moreover, the integrable majorants from Assumption (ref) together with the smoothness of $\psi$ allow to interchange the integral and the derivative of $\psi$.

assumpMAs $t\to\infty$, \begin{enumerate} • $M(t)\mathcal{I}^{-1}(t,\bm{\rho}_0)$ converges in probability to a positive semidefinite matrix; • $\int_0^{t}\left\{\mathcal{I}^{-1/2}(t,\bm{\rho}_0)\mathcal{K}(z,\bm{\rho}_0)\mathcal{I}^{-1/2}(t,\bm{\rho}_0)\right\}^2\mbox{d} z\to\bm{0}$. \end{enumerate}

To check whether Assumption (ref) (i) holds, one needs, for instance, to verify the following two relations: $\mathcal{I}^{-1}(t,\bm{\rho}_0)\mathsf{E} M(t)=\mathcal{I}^{-1}(t,\bm{\rho}_0)\Psi(t;\bm{\rho}_0)$ converges to a positive semidefinite matrix, which may also be a zero matrix; and $\mathsf{Var} \big\{\left(\mathcal{I}^{-1}(t,\bm{\rho}_0)\right)_{i,j}M(t)\big\}=\left(\mathcal{I}^{-1}(t,\bm{\rho}_0)\right)_{i,j}^2\Psi(t;\bm{\rho}_0)\to 0$ as $t\to\infty$ for all $i,j=1,\ldots,q$. By the Cauchy-Schwarz inequality and equations (ref)--(ref) from the proof of the consequent Theorem (ref), Assumption (ref) (ii) is satisfied if for all $j,k,\ell,m=1,\ldots,q$ holds

align*[align* omitted — 389 chars of source]

At first sight, technical Assumptions (ref) and (ref) have also a practical impact on admissibility of the intensity function $\psi$ and on the amount of information about the parameter $\bm{\rho}$ contained in the process $M$. Basically, the intensity $\psi$ has to be sufficiently smooth and adequately regular with respect to the information matrix function $\mathcal{I}$. In practice, we do not allow for too `wild' and too quickly changing behavior of the process of reporting times.

theorem[Consistency I] Under Assumptions (ref)--(ref), \[ \mathcal{I}^{1/2}(t,\bm{\rho}_0)\left( \widehat{\bm{\rho}}-\bm{\rho}_0\right)=-\mathcal{I}^{-1/2}(t,\bm{\rho}_0)\sum_{i=1}^{M(t)}\partial_{\bm{\rho}_0}h\left(Z_i;\bm{\rho}_0,t\right)+o_{\mathsf{P}}(1),\quad t\to\infty. \]

Let us discretize the `continuous' time $t\in\mathbb{R}_0^+$ for the process $\{M(t)\}_{t\geq 0}$ in a way that one observes $M$ only at all discrete time points $a\in\mathbb{N}$. This is indeed in concordance with the nature of our practical problem, where we evaluate the number of reported claims at the end of the calendar year represented by a discrete value of $a$.

Additionally, the next Lindeberg condition can extend the assertion of Theorem (ref).

assumpM$\lim_{a\to\infty}\sum_{i=1}^{a}\mathsf{E}\left({\boldsymbol d}^{\top}\boldsymbol Y_i\right)^2\mathbbm{1}\{|{\boldsymbol d}^{\top}\boldsymbol Y_i|\geq \varepsilon\|{\boldsymbol d}\|_2\}=0$ for all ${\boldsymbol d}\in\mathbb{R}^{q}$ and $\varepsilon>0$, where $\boldsymbol Y_i:=\mathcal{I}^{-1/2}(a,\bm{\rho}_0)\int_{i-1}^{i}\left\{\partial_{\bm{\rho}_0}\log\psi(z;\bm{\rho})\right\}\left(\mbox{d} M(z)-\psi(z;\bm{\rho}_0)\mbox{d} z\right)$.

For practical verification purposes, one can assume a version of the Lyapunov condition instead of the Lindeberg one. For instance, for all ${\boldsymbol d}\in\mathbb{R}^{q}$, there exists $\delta>0$ such that $\lim_{a\to\infty}\sum_{i=1}^{a}\mathsf{E}\left|{\boldsymbol d}^{\top}\boldsymbol Y_i\right|^{2+\delta}=0$. On one hand, the Lyapunov condition is more restrictive than the Lindeberg one, on the other hand, it is easier to verify.

corollary[Asymptotic normality I] Under Assumptions (ref)--(ref), \[ \mathcal{I}^{1/2}(a,\bm{\rho}_0)\left( \widehat{\bm{\rho}}-\bm{\rho}_0\right)\xrightarrow[a\to\infty]{\mathsf{D}}\mathsf{N}_{q}\left(\bm{0},{\boldsymbol I}\right). \]

Let us consider cases of the constant and exponential intensity function.

exampleIntensity $\psi(z;\bm{\rho})=\rho$, which corresponds to a homogeneous Poisson process. Here, ${\boldsymbol R}=(0,\infty)$ and the cumulative intensity is $\Psi(t,\rho)=\rho t$. The log-likelihood function is $\ell(\rho;\bm{Z},t)=M(t)\log\rho-\rho t$. The ML estimator is $\hat{\rho}=M(t)/t$, the information number becomes $\mathcal{I}(t;\rho_0)=t/\rho_0$, and $\mathcal{K}(t;\rho_0)=-\rho_0^{-3/2}$. Assumption (ref) is easily verifiable. For the Lyapunov condition, the choice of $\delta=1$ leads to $\sum_{i=1}^{a}\left(\frac{\rho_0}{a}\right)^{3/2}\mathsf{E}\left|\int_{i-1}^i\frac{1}{\rho_0}\mbox{d} M(z)-\int_{i-1}^i\mbox{d} z\right|^3=\sqrt{\frac{\rho_0^{3}}{a}}\mathsf{E}\left|\frac{Z_{i}-Z_{i-1}}{\rho_0}-1\right|^3=\sqrt{\frac{\rho_0^{3}}{a}}(12\mathrm{e}^{-1}-2)\to 0$ as $a\to\infty$, because the interarrival times $\{Z_i-Z_{i-1}\}_i$ have the exponential distribution with parameter $1/\rho_0$ (i.e., its expectation equals $\rho_0$). Hence, $\sqrt{a/\rho_0}\left( \widehat{\rho}-\rho_0\right)\xrightarrow[a\to\infty]{\mathsf{D}}\mathsf{N}\left(0,1\right)$.
exampleIntensity $\psi(z;\bm{\rho})=\exp\{\rho_1+\rho_2z\}$. Here, ${\boldsymbol R}=(0,\infty)\times (0,\infty)$. The log-likelihood function is $\ell(\bm{\rho};\bm{Z},t)=\rho_1M(t)+\rho_2\sum_{i=1}^{M(t)}Z_i-\mathrm{e}^{\rho_1}\left(\mathrm{e}^{\rho_2 t}-1\right)/\rho_2$. The ML estimator of the parameter $\rho_2$ can be obtained as a solution of $\sum_{i=1}^{M(t)}Z_i+M(t)/\hat{\rho}_2-tM(t)/\left(1-\mathrm{e}^{-\hat{\rho}_2t}\right)=0$ and the ML estimator of the parameter $\rho_1$ comes from $\hat{\rho}_1=\log\left\{\hat{\rho}_2M(t)/\left(\mathrm{e}^{\hat{\rho}_2t}-1\right)\right\}$. Consequently, \[ \mathcal{I}(t;\bm{\rho})=\left(\begin{matrix} \mathrm{e}^{\rho_1} \left(\mathrm{e}^{\rho_2 t}-1\right)/\rho_2 & \mathrm{e}^{\rho_1} \left\{\mathrm{e}^{\rho_2 t} \left(\rho_2 t-1\right)+1\right\}/\rho_2^2\\ \mathrm{e}^{\rho_1} \left\{\mathrm{e}^{\rho_2 t} \left(\rho_2 t-1\right)+1\right\}/\rho_2^2 & \mathrm{e}^{\rho_1} \left[\mathrm{e}^{\rho_2 t} \left\{\rho_2 t \left(\rho_2 t-2\right)+2\right\}-2\right]/\rho_2^3 \end{matrix}\right), \] which can be easily proved to be positive definite for all $t>0$ and any $\bm{\rho}\in{\boldsymbol R}$. Assumption (ref) together with the Lyapunov condition can be checked as well. Hence, $\mathcal{I}^{1/2}(a;\bm{\rho}_0)\left(\widehat{\bm{\rho}}-\bm{\rho}_0\right)\xrightarrow[a\to\infty]{\mathsf{D}}\mathsf{N}_2\left(\bm{0},{\boldsymbol I}\right)$.

An example of the intensity function directly used in the consequent practical analysis of our data for modeling the reporting times of bodily injury claims, see also Figure (ref) (left panel, green line), is given below.

exampleIntensity $\psi(z;\bm{\rho})=\exp\{\rho_1+\rho_2\log z+\rho_3\cos\left(2\pi z/\rho_5\right)+\rho_4\sin\left(2\pi z/\rho_5\right)\}$. The above formulated assumptions are satisfied for a particular open convex ${\boldsymbol R}\subseteq\mathbb{R}^5$. The defined entities are not presented here due to their voluminous forms.
figure[figure omitted — 401 chars of source]

Besides that, the next example is used in the data analysis for the reporting times of material damage claims, cf. Figure (ref) (right panel, green line).

exampleIntensity $\psi(z;\bm{\rho})=\exp\{\rho_1+\rho_2 z+\rho_3 z^2+\rho_4\cos\left(2\pi z/\rho_6\right)+\rho_5\sin\left(2\pi z/\rho_6\right)\}$. The required assumptions are again satisfied, but the above defined entities are not presented here due to their complicated and voluminous forms.

Suitability of Examples (ref) and (ref) for the practical analysis is illustrated in Figure (ref), where the observed and fitted cumulative intensities (corresponding to the theoretical cumulative intensity $\Psi(t;\bm{\rho})$) are compared. The deviations between them are minor. Let us recall that our estimation of the reporting dates' intensity is based on the data up to the end of year 2015. Extrapolation of the estimated cumulative intensity for the `future' year 2016 also nicely mimics the known reality from year 2016 (not used for estimation). The estimated underlying intensity is depicted as well.

Finally, it would be natural to characterize the number of new arriving claims with respect to the intensity of the process $M$.

proposition[Infinite number of renewals] If Assumption (ref) holds and $\lim_{t\to\infty}\Psi(t;\bm{\rho})=\infty$, then $\mathsf{P}\left\{\lim_{t\to\infty}M(t)=\infty\right\}=1$.

This proposition reveals that the divergent cumulative intensity $\Psi(t;\bm{\rho})$ assures that there are still new claims being reported with probability one.

Number of payments

Let us recall that the number of payments corresponding to the $i$th claim till time point $t$ is denoted by $N_i(t)$. So, we possess panels of count processes $\{N_i(t)\}_{t>0}$ for $i=1,\ldots,M(t)$ that can be represented as $\{\bm{N}(t)\}_{t>0}$. Thus, the claim notifications together with the claim payments can be viewed as a marked Poisson process with Poisson processes as marks \[ \{\{M(t)\}_{t\geq 0},\{\bm{N}(t)\}_{t\geq 0}\}. \] In practice, we observe the counting process $\{N_i(t)\}_{t>0}$ through $\{U_{i,1},\ldots,U_{i,N_i(t)}\}$ for $i=1,\ldots,M(t)$, where the number of payments for the $i$th claim is denoted by $N_i(t)$ and $U_{i,k}$ is the time of the $k$th payment within the $i$th claim for $k=1,\ldots,N_i(t)$. Moreover, the amount of the $k$th payment for the $i$th claim paid at time $U_{i,k}$ is represented by $X_{i,k}$, which is going to be modeled in Subsection (ref).

assumpNThe ordered payment times $\{U_{i,1},U_{i,2},\ldots\}$ of the $i$th loss are arrival times of a non-homogeneous Poisson process $\{N_i(t)\}_{t\geq 0}$. Processes $\{N_i(t)\}_{t\geq 0}$, $i=1,2,\ldots$ are independent having parametric intensities $\lambda(t,Z_i;\bm{\theta})$ such that $N_i(t)=\sum_{k=1}^{\infty}\mathbbm{1}\{U_{i,k}\leq t\}$, $\bm{\theta}\in{\boldsymbol P}\subseteq\mathbb{R}^p$, and ${\boldsymbol P}$ is an open convex set.

Since $U_{i,k}\geq Z_i$ (i.e., the payment times come after the reporting time), the corresponding density $\lambda$ has to be constant zero up to the reporting date $Z_i$. Alternatively, one can think of a `restarted' process $\tilde{N}_{Z_i}(\tau)=\sum_{k=1}^{\infty}\mathbbm{1}\{U_{i,k}-Z_i\leq\tau\}$ with an `internal' time $\tau$ of the claim $i$ after its reporting time $Z_i$.

Note that the processes $\{N_i(t)\}_{t\geq 0}$, $i=1,2,\ldots$ do not have to be identically distributed, because of possible different effect of the reporting date. Here, the intensities $\lambda(t,Z_i;\bm{\theta})$ can be considered as payment frequencies. However, the common parameter $\bm{\theta}$ is assumed to be shared by different intensity functions $\lambda$'s. In contrast to Fawless1987, we do not assume a specific product form of the intensity $\lambda$ and, moreover, we allow the processes $N_i$ to depend on the process $M$ via the reporting times $Z_i$'s as the marks' locations. To the best of our knowledge, we are not aware of any previous work dealing with inference for the marked non-homogeneous Poisson process with marks being non-homogeneous Poisson processes.

In order to estimate the unknown parameter $\bm{\theta}$, we again use the ML approach for the arrival times, which can be considered as an extension of the case for a single realization of the Poisson process. Such a framework can be extended for several independent non-homogeneous Poisson processes, where the likelihood is as follows

equation*[equation* omitted — 225 chars of source]

where $\bm{Z}=(Z_1,\ldots,Z_{M(t)})^{\top}$ can be also viewed as the covariates (regressors) of the intensity $\lambda$ having corresponding realizations $\bm{z}$. One should bear in mind that the whole information about the process $\{M(t)\}_{t\geq 0}$ is included in the sequence $\{Z_1,Z_2,\ldots\}$. Furthermore, the intensity function may be decomposed, for instance, as

equation[equation omitted — 92 chars of source]

such that $\bm{\theta}=(\bm{\nu}^{\top},\bm{\eta}^{\top})^{\top}$, $\lambda_0(\tau;\bm{\nu})$ is a baseline intensity function, where $\lambda_0(\tau;\bm{\nu})=0$ for $\tau<0$, and $f(Z_i;\bm{\eta})$ is a parametric covariate function introducing the effects of the covariates $Z_i$'s.

To obtain the ML estimator of $\bm{\theta}$, one has to maximize the log-likelihood function

equation[equation omitted — 212 chars of source]

The true value $\bm{\theta}_0\in{\boldsymbol P}$ of the unknown vector parameter $\bm{\theta}$ is supposed to uniquely maximize $\mathsf{E}_{\bm{\theta}}\ell\{\bm{\theta};\bm{N}(t),M(t)\}$. For $\bm{U}_i:=(U_{i,1},\ldots,U_{i,N_i(t)})^{\top}$, let us define \[ g_i(\bm{U}_i;\bm{\theta},t):=\int_{Z_i}^{t}\lambda(\tau,Z_i;\bm{\theta})\mbox{d} \tau-\sum_{k=1}^{N_i(t)}\log\lambda(U_{i,k},Z_i;\bm{\theta}), \] which means that

equation[equation omitted — 142 chars of source]

Assume the following to hold with respect to the $g_i(\bm{U}_i;\bm{\theta},t)$ functions in order to obtain consistent and asymptotically normal estimators. Next assumption, being analogous to Assumption (ref), secures the existence of the unique solution.

assumpN$g_i(\bm{u};\bm{\theta},t)$ are convex in $\bm{\theta}\in{\boldsymbol P}$ for all $0<u_1<\ldots<u_n<t$ and $n\in\mathbb{N}$.

At a very first sight, further $\mathscr{N}$-assumptions might be considered as the copied $\mathscr{M}$-assump\-tions mutatis mutandis. There is, however, an additional layer of randomness present for the marks $N_i$'s. For instance, deterministic integrals become stochastic ones. Furthermore, the assumptions regarding the process $M$ are not just replaced by the assumptions for the processes $N_i$'s. There are indeed several additional assumptions regarding the marks $N_i$'s added to some assumptions for the original underlying process $M$. Therefore, one can neither simplify nor unify the $\mathscr{M}$- and $\mathscr{N}$-assumptions. Next assumption, being similar to Assumption (ref), controls the almost sure boundedness of the derivatives of the intensity functions for all $Z_i$.

assumpN$\frac{\partial^2}{\partial\bm{\theta}\partial\bm{\theta}^{\top}}\lambda(\cdot,Z_i;\bm{\theta}):\,\mathbb{R}^+\to\mathbb{R}^{p\times p}$ are continuous for all $\bm{\theta}\in{\boldsymbol P}$ and there exist Lebesgue-integrable functions $m_{1,i,j}$ and $m_{2,i,j,k}$ such that \[ \left|\frac{\partial}{\partial\theta_j}\lambda(t,Z_i;\bm{\theta})\right|\leq m_{1,i,j}(t)\quad\mbox{and}\quad\left|\frac{\partial^2}{\partial\theta_j\partial\theta_k}\lambda(t,Z_i;\bm{\theta})\right|\leq m_{2,i,j,k}(t) \] almost surely, for all $\bm{\theta}\in{\boldsymbol P}$, $i\in\mathbb{N}$, almost every $t>0$, and $j,k=1,\ldots,p$.

Let us define a cumulative intensity $\Lambda(t,Z_i;\bm{\theta}):=\int_{Z_i}^t\lambda(\tau,Z_i;\bm{\theta})\mbox{d}\tau$, an information matrix $\mathcal{J}(t;\bm{\theta}_0):=\mathsf{E}\sum_{i=1}^{M(t)}\mathcal{J}_i(t;\bm{\theta}_0)$, where $\mathcal{J}_i(t;\bm{\theta}_0):=\int_{Z_i}^t\frac{\{\partial_{\bm{\theta}_0}\lambda(\tau,Z_i;\bm{\theta})\}^{\otimes 2}}{\lambda(\tau,Z_i;\bm{\theta}_0)}\mbox{d} \tau$, and a matrix $\mathcal{L}_i(t,\bm{\theta}_0):=\frac{1}{\sqrt{\lambda(t,Z_i;\bm{\theta}_0)}}\left[\partial^2_{\bm{\theta}_0}\lambda(t,Z_i;\bm{\theta})-\frac{\{\partial_{\bm{\theta}_0}\lambda(t,Z_i;\bm{\theta})\}^{\otimes 2}}{\lambda(t,Z_i;\bm{\theta}_0)}\right]$ for $t>0$. The forthcoming assumption differs from the analogous one (ref) via averaging over all the payments $Z_i$.

assumpNAs $t\to\infty$, \begin{enumerate} • $M(t)\mathcal{J}^{-1}(t,\bm{\theta}_0)$ converges in probability to a positive semidefinite matrix; • $\mathsf{E}\sum_{i=1}^{M(t)}\int_{Z_i}^{t}\left\{\mathcal{J}^{-1/2}(t,\bm{\theta}_0)\mathcal{L}_i(\tau,\bm{\theta}_0)\mathcal{J}^{-1/2}(t,\bm{\theta}_0)\right\}^2\mbox{d} \tau\to\bm{0}$. \end{enumerate}

Analogous discussions like after Assumptions (ref)--(ref) might be carried out regarding Assumptions (ref)--(ref). Briefly and informally, the intensity functions $\lambda$'s representing the behavior of the payment times' processes $N_i$'s are supposed to be sufficiently smooth and adequately regular with respect to the amount of information about the parameter $\bm{\theta}$ contained in $N_i$'s.

theorem[Consistency II] Under Assumptions (ref) and (ref)--(ref), \[ \mathcal{J}^{1/2}(t,\bm{\theta}_0)\left(\widehat{\bm{\theta}}-\bm{\theta}_0\right)=-\mathcal{J}^{-1/2}(t,\bm{\theta}_0)\sum_{i=1}^{M(t)}\partial_{\bm{\theta}_0}g_i\left(\bm{U}_i;\bm{\theta}_0,t\right)+o_{\mathsf{P}}(1),\quad t\to\infty. \]

Again, let us discretize the `continuous' time $t\in\mathbb{R}_0^+$ for the processes $\{N_i(t)\}_{t\geq 0}$ in a way that one observes $N_i$ only at discrete time points $a\in\mathbb{N}$, e.g., status at the closed calendar years. The next Lindeberg condition is used to extend the assertion of Theorem (ref) in order to derive asymptotic normality of $\widehat{\bm{\theta}}$.

assumpN$\lim_{a\to\infty}\sum_{i=1}^{a}\mathsf{E}\left({\boldsymbol d}^{\top}\mathcal{Y}_i\right)^2\mathbbm{1}\{|{\boldsymbol d}^{\top}\mathcal{Y}_i|\geq \varepsilon\|{\boldsymbol d}\|_2\}=0$ for all ${\boldsymbol d}\in\mathbb{R}^{p}$ and $\varepsilon>0$, where $\mathcal{Y}_i:=\mathcal{J}^{-1/2}(a,\bm{\theta}_0)\int_{j-1}^{j}\int_{z}^{a}\left\{\partial_{\bm{\theta}_0}\log\lambda(\tau,z;\bm{\theta})\right\}(\mbox{d}\tilde{N}_z(\tau-z)-\lambda(\tau,z;\bm{\theta}_0)\mbox{d}\tau)\mbox{d} M(z)$.

The Lyapunov condition can be assumed as well. For instance, for all ${\boldsymbol d}\in\mathbb{R}^{p}$, there exists $\delta>0$ such that $\lim_{a\to\infty}\sum_{i=1}^{a}\mathsf{E}\left|{\boldsymbol d}^{\top}\mathcal{Y}_i\right|^{2+\delta}=0$.

corollary[Asymptotic normality II] Under Assumptions (ref) and (ref)--(ref), \[ \mathcal{J}^{1/2}(a,\bm{\theta}_0)\left( \widehat{\bm{\theta}}-\bm{\theta}_0\right)\xrightarrow[a\to\infty]{\mathsf{D}}\mathsf{N}_{p}\left(\bm{0},{\boldsymbol I}\right). \]

The simplest situation is that each $\{N_i(t)\}_{t\geq 0}$ is a homogeneous Poisson process after the reporting time $Z_i$ having a constant common intensity $\theta>0$ for all $i$'s.

exampleIntensity $\lambda(\tau,Z_i;\theta)=\theta\mathbbm{1}\{\tau\geq Z_i\}$ with ${\boldsymbol P}=(0,\infty)$. The log-likelihood function is $\ell\{\theta;\bm{N}(t),M(t)\}=(\log\theta)\sum_{i=1}^{M(t)} N_i(t)-\theta\sum_{i=1}^{M(t)}(t-Z_i)$. The ML estimator becomes $\hat{\theta}=\frac{\sum_{i=1}^{M(t)}N_i(t)}{tM(t)-\sum_{i=1}^{M(t)}Z_i}$ and the information number is $\mathcal{J}(t;\theta_0)=\left\{t\Psi(t;\bm{\rho}_0)-\int_0^tz\psi(z;\bm{\rho}_0)\mbox{d} z\right\}/\theta_0$. Assumption (ref) together with the Lyapunov condition can be checked as well. Hence, \[ \sqrt{\frac{a\Psi(a;\bm{\rho}_0)-\int_0^az\psi(z;\bm{\rho}_0)\mbox{d} z}{\theta_0}}\left( \widehat{\theta}-\theta_0\right)\xrightarrow[a\to\infty]{\mathsf{D}}\mathsf{N}\left(0,1\right). \] Moreover, if the underlying intensity of the Poisson process $M$ is $\psi(t;\rho)=\rho$ for $t\geq 0$ (i.e., homogeneous Poisson process), then $\sqrt{\frac{\rho_0}{2\theta_0}}a\left( \widehat{\theta}-\theta_0\right)\xrightarrow[a\to\infty]{\mathsf{D}}\mathsf{N}\left(0,1\right)$.

Another example contains a baseline intensity, which is motivated by CKSS1991.

exampleIntensity $\lambda(\tau,Z_i;\bm{\nu},\eta)=\nu_1\nu_2(\tau-Z_i)^{\nu_1-1}\exp\{\eta Z_i\}\mathbbm{1}\{\tau\geq Z_i\}$. Here, ${\boldsymbol P}=(0,\infty)^2\times\mathbb{R}$. The log-likelihood becomes \begin{multline*} \ell\{\bm{\nu},\eta;\bm{N}(t),M(t)\}\\=\sum_{i=1}^{M(t)}\left\{N_i(t)\log(\nu_1\nu_2)+(\nu_1-1)\sum_{k=1}^{N_i(t)}\log(U_{i,k}-Z_i)+\eta Z_iN_i(t)-\nu_2\left(t-Z_i\right)^{\nu_1}\exp\{\eta Z_i\} \right\}. \end{multline*} The ML estimator has to be computed numerically. The above formulated assumptions are satisfied for the particular open convex ${\boldsymbol P}\subseteq\mathbb{R}^3$. The defined entities are not presented here due to their complicated forms.

The next example of the intensity function is directly used in the consequent practical analysis of our data for modeling the payment times of bodily injury as well as material damage claims.

example$\lambda(\tau,Z_i;\bm{\nu},\bm{\eta})=\exp\{\nu_1+\nu_2(\tau-Z_i)+\eta_1\cos\left(2\pi Z_i/\eta_3\right)+\eta_2\sin\left(2\pi Z_i/\eta_3\right)\}$. The defined entities are again not presented here due to their voluminous forms.

Reporting delay

Since the reporting delays correspond to different claims from different accidents, independence between the reporting delays $W_i$'s for different contracts is assumed. However, the distribution of the reporting delays is allowed to change with respect to the reporting time $Z_i$. The reporting delays seem to become shorter and shorter, which can be explained by a possibility to report an accident over the internet and even by a denser net of the insurance company branches. So, given $Z_i$, $W_i$ has a parametric conditional density $f_{W}\{\cdot;w(Z_i,\bm{\vartheta})\}$, where $\bm{\vartheta}\in\mathbb{R}^r$. Note that $\{W_i\}_{i\in\mathbb{N}}$ are not identically distributed, which allows, for instance, to assume time-varying distributions for the $W_i$'s through the function $w(\cdot,\bm{\vartheta})$. Using similar arguments as in HjortPollard2011, one has consistency and asymptotic normality for the ML estimator $\widehat{\bm{\vartheta}}$.

A variety of parametric distributions are suitable for the analysis, particularly those with the shapes similar to the ones provided by log-normal, Weibull, or Gamma distributions. All of them have similar performance, whereas the importance lies in the time-varying parameters. Although in the rest of the study, we concentrate ourselves purely on the log-normal distribution, let us briefly recall all three mentioned densities

align[align omitted — 383 chars of source]

For the Gamma and Weibull distributions, parameters $c$ and $d$ are called the shape and scale, respectively, for the log-normal distribution $c$ and $d$ are the mean and standard deviation of the distribution on the log scale, respectively.

Taking into account the dependency between the accident date $Z_i$ and the reporting delay $W_i$ and, additionally to that, allowing for possible seasonal behavior, we consider truncated Fourier series for the parameters of the conditional distributions with

align*[align* omitted — 274 chars of source]

in the form of

align[align omitted — 564 chars of source]

where $52$ is the number of weeks in one year and 7 is the number of days in one week. Later on for considering models with increasing flexibility, we discuss constant models with $L = 0$ and $\beta_c=\beta_d=0$, linear models with $L=0$, and models with one ($L = 1$) or two ($L= 2$) seasonality patterns. The ML estimator of $\bm{\vartheta} = (\bm{\vartheta}_1^\top, \bm{\vartheta}_2^\top)^\top$ is obtained as

align[align omitted — 317 chars of source]

The condition $c(Z_i,\bm{\vartheta}_1)>0$ is considered only in case when the assumed distribution $f_W$ is Weibull or Gamma. As initial values in the iterative maximization procedure, we took the parameter values that at best fit the piecewise constant weekly averaged values in the least squares sense.

Figure (ref) shows the estimation results for $f_W$ being the log-normal distribution.

figure[figure omitted — 535 chars of source]

The left panel corresponds to the bodily injury claims, where the right one to the material damage claims, on the bottom panel we show the number of accidents as the function of the reporting date. In the two upper panels, we depict the location parameters $c(z, \widehat{\bm{\vartheta}}_1)$ and two middle panels scale parameters $d(z, \widehat{\bm{\vartheta}}_2)$ over different weeks of the reporting date. Different colors represent different complexities of the models used: for cyan we have a constant unconditional model with $L = 0$ and $\beta_c=\beta_d=0$; pink uses only linear temporal dependency with $L = 0$; green and blue lines have one ($L = 1$) and two ($L = 2$) levels of seasonality. The most flexible density with $L = 2$ has been used in the final omnibus model. With red we depict values extrapolated to the data not used in the estimation, namely years 2003 and 2016. This analysis shows strong temporal dependence of the parameters on the distribution of $W_i$.

In order to validate our results from the fitted model, Figure (ref) compares the observed and predicted quarterly averaged reporting (waiting) delays in days for the bodily injury claims as well as for the material damage claims.

figure[figure omitted — 395 chars of source]

Payment amounts

Denote the $j$th payment amount for the $i$th claim by $X_{i,j}$, where $j=1,\ldots,N_i(t)$. The $X_{i,j}$'s are independent over all $j$'s as well as all $i$'s. This independency assumption can be easily relaxed and any other times series model (e.g., an autoregressive model) can be used instead. However, our empirical findings imply independency. Given $Z_i$, $X_{i,j}$ has a parametric conditional density $f_{X}\{\cdot;v(Z_i,{\boldsymbol\varsigma})\}$, where ${\boldsymbol\varsigma}\in\mathbb{R}^s$ and the function $v(\cdot,{\boldsymbol\varsigma})$ introduces the time-varying effects of $Z_i$'s. Moreover, the reporting delays $W_i$'s are also supposed to be independent from the payment amounts $X_{i,j}$'s, as well as all the payment amounts from the same claim are independent among each other. These assumptions are based on the preliminary empirical analysis of the pairwise relationships between the waiting time ($W_i$) and the first ($X_{i, 1}$), second ($X_{i, 2}$), third ($X_{i, 3}$), and fourth ($X_{i, 4}$) claim payment amounts shown in Figure (ref). For this plot, the data are transformed by the estimated cumulative distribution functions (cdf) $\hat F_X$ and $\hat F_W$ obtained from the plugged-in densities $\hat f_X(\cdot)\equiv f_X\{\cdot,v(Z_i,\widehat{\boldsymbol{\varsigma}})\}$ and $\hat f_W(\cdot)\equiv f_W\{\cdot,w(Z_i,\widehat{\boldsymbol{\vartheta}})\}$, respectively. Further, the data are transformed via the quantile function of the standard normal distribution. If the distributional assumptions are correct and the payments and delays are independent, the bivariate kernel density estimates and scatterplots of the transformed data should be suggestive of circular shapes as clearly visible in Figure (ref).

figure[figure omitted — 806 chars of source]

Using similar arguments as in HjortPollard2011, one can prove consistency and asymptotic normality for the ML estimator $\widehat{\boldsymbol\varsigma}$. The procedure for modeling the claim payments closely resembles the procedure mentioned when modeling the reporting delays in Subsection (ref). For the modeling of the payment amounts, we also considered more complex models, where previous payments were included as exogenous variables. This however did not bring any improvements.

Figure (ref) presents the time-varying parameters of log-normal distribution of bodily injury as well as material damage claims for the first payment in yellow (fully flexible model with two seasonal periods) and in grey (separately for each week). Other curves present parameters for all the payment amounts pooled together. As there is no much difference and results for the pooled models seemed to be more stable, we concentrate in the later only on the pooled ones, namely the most flexible with two periods and a linear trend.

figure[figure omitted — 591 chars of source]

Practical application and empirical results

To numerically illustrate the performance of our method, we use two data sets---bodily injury and material damage claims (cf. motivation and data description in Section (ref)). Let us recall that data from the last available year 2016 are used only for back-testing and comparison with the predicted results. Furthermore, our `micro' (granular, claim-by-claim) approach is also compared with a traditional standard actuarial technique---bootstrap chain-ladder EV1999---in combination with linear extrapolation of the reported claims in the next year. This `macro' approach is based on aggregation of data and, hence, it disregards the information about the policy and the claim's development. We refer to it from now on as the aggregated method.

Firstly, the previously described estimation procedures (Section (ref)) provide parameter estimates of our omnibus model. Secondly, a Monte Carlo prediction technique is involved in order to generate (simulate) the future claims' developments. In essence, prediction for a distribution of the total payments in year 2016, for which we also possess the real paid claim amounts.

Parameter estimation

All the estimates are obtained through the ML approach, which guaranties a proper stochastic inference. For the case of densities, it is widely known that under some regularity conditions the ML estimators are consistent and asymptotically normal. For the case of intensities, we have proved consistency and asymptotic normality of these ML estimators. Consequently, one can plug-in the estimated parameters into the parametric forms of the densities and intensities present in our micro model in order to have predicted (fitted) intensities of the reporting dates/payment dates and densities of the reporting delays/payment amounts. They are going to be used for the simulation of the future payments (dates and amounts).

In particular, the intensity function for modeling the reporting times of bodily injury claims comes from Example (ref) and in case of material damage claims from Example (ref). The reporting delay and the claim payments are modeled as in relations (ref)--(ref). And the intensity functions for the payment times of bodily injury as well as material damage claims come from Example (ref).

Monte Carlo predictions

Basically for each claim being reported up to the future time point $b$ from Figure (ref) (e.g., end of the next calendar year), one needs to simulate payment dates and corresponding payment amounts, which are going to be summed in each simulation's run. These sums of payments give us the simulated (empirical) predictive distribution of the total future payments. Hence, for the next year (in a general future time window $(a,b]$), we need to simulate the new payment dates for all already reported claims as well as the payment dates corresponding to the incurred but not reported claims. Consequently, it is requisite to generate a corresponding payment amount for every payment time within the time interval $(a,b]$.

Let us realize that we need to generate many realizations of the non-homogeneous Poisson process for each Monte Carlo simulation's run. To simulate the non-homoge\-neous Poisson process, we suggest to rely on the thinning algorithm by LS1979. The main reason for choosing this way to generate enormous number of realizations of the non-homogeneous Poisson process is that this approach can be applied to any rate function without the necessity of numerical integration or simulation of Poisson variables.

Primary goal: Prediction of distribution of the future cash flows

Our primary target is to predict the distribution of the total payment amounts within the future time period. Such a prediction is going to be obtained through Procedure (ref).

algorithm[algorithm omitted — 2,662 chars of source]

The predicted distribution of the forthcoming payment amounts within the next year is graphically displayed in Figure (ref). Here, our micro approach is compared with the traditional method based on data aggregation. Moreover, the true cumulative amount of the one-year ahead payments is depicted in order to judge the point prediction's precision.

figure[figure omitted — 705 chars of source]

There are two general and, from a practical point of view, very important findings with respect to the prediction of the future total payment amounts. First, our claim-by-claim based method is more precise in point prediction to the `unknown' true value compared to the traditional technique based on aggregated data. And this holds for both lines of business. Second, our micro approach provides less volatile predicted distribution, e.g., in terms of the coefficient of variation.

Secondary goal: Back-prediction of the truncated occurrence times

Our secondary practical target is to back-predict the accident dates of the claims, which are truncated due to the reporting delay. We are indeed not aware of so-called incurred but not reported claims and the insurance company needs to back-predict these claims, which have already occurred, but have not been reported yet. This can be reached via Procedure (ref).

algorithm[algorithm omitted — 983 chars of source]

The counts of the back-predicted accident dates for the latest year as well as the counts of the predicted accident dates for the next year are visualized in Figure (ref).

figure[figure omitted — 600 chars of source]

It is of utmost importance for the insurance company to have information regarding the accidents, which have not been reported yet. Figure (ref) reveals triplets of bars, that nicely accommodate the problem of truncated data in our setup. The middle horizontal bar in each triplet corresponds to the known number of accident dates based on the database up to year 2015. The left bar in the triplet stands again for the known number of accident dates, although, coming from the database up to year 2016. Therefore, the left bar has to be higher than the middle one, because additional claims can be reported within the calendar year 2016 (and they can occur in 2016 or even in previous years). The right bar in the triplet represents our prediction. It is supposed to be slightly higher even than the left bar, because there can occur additional accidents even before the end of year 2016 that are not going to be reported till the end of year 2016.

Conclusions

Micro forecasting is in general a stochastic prediction method for future losses/costs relying on the individual developments of the recorded historical events. Our prediction approach is capable to model the probabilistic behavior of the future losses' occurrences, the occurrences of the incurred but not reported losses, the lengths of the reporting delays, and the frequency and severity of the loss payments in time. This is indeed sufficient for prediction of the future cash-flows for a predetermined time horizon.

We employ the micro prediction technique in claims reserving. To meet all future insurance claims rising from policies, it is requisite to quantify the outstanding loss liabilities. Here, utility for solvency of the insurance company is developed. And, clearly, valuation of the reserving risk in insurance is not the only area of empirical economics, where the proposed methodology can be applied as documented by several case examples.

Quantifying reserving risk in non-life insurance inadvertently yields to a theoretical framework of the marked non-homogeneous Poisson process with non-homogeneous Poisson processes as marks. It can be viewed as an infinitely stochastic Poisson process and, consequently, a proper statistical inference relying on simple and verifiable assumptions is derived.

Acknowledgements

The research of Mat\'{u}\v{s} Maciak and Michal Pe\v{s}ta was supported by the Czech Science Foundation project GA\v{C}R No. 18-01781Y.