Jilling F Bergmann

Joep De Bruijne

Daphne M Hotho

Robert J De Knegt

Andre Boonstra

Christine J Weegink

Andre A Van Vliet

Jeroen Van De Wetering

Simon P Fletcher

Lisa A Bauman

Jilling ; Bergmann

Lisa ; Bauman

Pharmaceuticals, Inc Anadys

Rahimy Mohamad ; Anadys

MD Harry L A Janssen
email: h.janssen@erasmusmc.nl

J F Bergmann

L A Bauman

M Rahimy

J R Appleman

are employees of J L Freddo

H W Reesink

Roche S P Fletcher

L A Bauman

Anadys J F Pharmaceuticals

J Bergmann

D M De Bruijne

A Hotho

C J Boonstra

J Weegink

A A Van De Wetering

Antiviral Activity of ANA773, an Oral Inducer of Endogenous Interferons that Acts Via TLR7, in Chronic Hepatitis C Patients

Keywords: toll-like receptor agonist, safety, pharmacokinetics, pharmacodynamics, immunological response TLR, toll-like receptor, IFN-α, interferon-α, HCV, hepatitis C virus, RNA, ribonucleic acid, RBV, ribavirin, SVR, sustained virological response, HIV, human immunodeficiency virus, HBV, hepatitis B virus, LLQ, lower limit of quantification, ELISA, enzyme-linked immunosorbent assay, ISG15, interferon-stimulated gene 15, OAS, 2',5'-oligoadenylate synthetase, RIA, radioimmunoassay, ANOVA,

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Chronic hepatitis C virus (HCV) infection is a major cause of liver cirrhosis and hepatocellular carcinoma. HCV-related end-stage liver disease is now the main indication for liver transplantation in North America and Western Europe [START_REF] Terrault | Treating hepatitis C infection in liver transplant recipients[END_REF]. Estimates suggest that there are 170 million HCV infected patients worldwide and 3 to 4 million people are newly infected each year [START_REF] Lauer | Hepatitis C virus infection[END_REF].

Approximately 80% of patients who become infected with HCV develop chronic hepatitis C (3).

Combination therapy with pegylated interferon-alfa (PEG-IFN-α) and ribavirin (RBV) results in sustained virologic response (SVR) in only approximately 40-50% of patients with genotype 1 [START_REF] Fried | Peginterferon alfa-2a plus ribavirin for chronic hepatitis C virus infection[END_REF][START_REF] Hadziyannis | Peginterferon-alpha2a and ribavirin combination therapy in chronic hepatitis C: a randomized study of treatment duration and ribavirin dose[END_REF][START_REF] Manns | Peginterferon alfa-2b plus ribavirin compared with interferon alfa-2b plus ribavirin for initial treatment of chronic hepatitis C: a randomised trial[END_REF]. Thus, there remains an urgent unmet medical need to offer new therapies that may eradicate HCV infection [START_REF] Bergmann | What is on the horizon for treatment of chronic hepatitis C?[END_REF]. The most likely agents that will enter clinical practice soon are the protease inhibitors telaprevir and boceprevir. Phase III clinical trials of these drugs combined with peginterferon and ribavirin in genotype 1 patients have shown SVR rates up to 80% in treatmentnaive patients and up to 65% in treatment-experienced patients [START_REF] Bacon | Boceprevir for previously treated chronic HCV genotype 1 infection[END_REF][START_REF] Jacobson | Telaprevir in combination with peginterferon and ribavirin in genotype 1 HCV treatment-naïve patients: final results of Phase 3 ADVANCE study[END_REF][START_REF] Sherman | Telaprevir in combination with peginterferon alfa2a and ribavirin for 24 or 48 weeks in treatment-naive genotype 1 HCV patients who achieved Antiviral activity of ANA773 in chronic HCV patients Final version 22JAN2011 an extended rapid viral response: final results of Phase 3 ILLUMINATE study[END_REF][START_REF] Poordad | Boceprevir for untreated chronic HCV genotype 1 infection[END_REF].

Immune responses against HCV are considered to play an important role in controlling and clearing chronic viral HCV infections [START_REF] Ghany | Diagnosis, management, and treatment of hepatitis C: an update[END_REF][START_REF] Boonstra | Immunology of hepatitis B and hepatitis C virus infections[END_REF]. Toll-Like Receptors (TLRs) are a family of pathogen-recognition receptors that activate the innate immune response [START_REF] Iwasaki | Toll-like receptor control of the adaptive immune responses[END_REF]. Stimulation of TLRs either directly or indirectly leads to (i) the release of multiple cytokines, including type I and type II IFNs, (ii) the induction of pathways and enzymes that destroy intracellular pathogens and (iii) the stimulation of antigen-presenting cells which may result in the activation of adaptive immunity. To date, ten functional TLRs have been identified in humans. Among these, TLR7 is of particular interest because small molecule agonists for this receptor have been identified [START_REF] Fletcher | Masked oral prodrugs of toll-like receptor 7 agonists: a new approach for the treatment of infectious disease[END_REF], and like TLR9 agonists, induces type III IFNs. Furthermore, there is precedent for this target in viral diseases. In a clinical proof-of-concept study, a TLR7 agonist, isatoribine, was shown to have an anti-viral effect in the treatment of chronic HCV infection [START_REF] Horsmans | Isatoribine, an agonist of TLR7, reduces plasma virus concentration in chronic hepatitis C infection[END_REF]. In addition, topical imiquimod is approved for the treatment of genital papillomavirus infections [START_REF] Hemmi | Small antiviral compounds activate immune cells via the TLR7 MyD88-dependent signaling pathway[END_REF]. Unfortunately, imiquimod is poorly tolerated when administered orally [START_REF] Savage | A phase I clinical trial of imiquimod, an oral interferon inducer, administered daily[END_REF] and is both extensively and rapidly metabolized when administered systemically, limiting systemic exposure of the drug.

ANA773 is an oral prodrug of a TLR7 agonist, developed for the treatment of patients with chronic HCV infection. The active metabolite of ANA773 induces multiple subtypes of endogenous IFN-α, induces IFN-dependent anti-HCV activity and functionally activates natural killer (NK) cells in vitro. Oral and intravenous administration of ANA773 induced systemic IFN-α production in Cynomolgus monkeys, with IFN-dependent responses being observed in both the periphery and the liver. The alternate-day dosing schedule in these studies were shown to induce desirable levels of T-cell activation and proliferation, and has been selected as the dosing schedule for clinical evaluation. A single and multiple dose escalation study (200-1600 mg/day every-other-day up to 4 administrations) of oral ANA773 in 41 healthy volunteers showed a dose-related increase of IFN-α and markers of IFN response. The occurrence of treatment related adverse events in this study were seen predominantly in the higher dose groups (1200 and 1600 mg/day) and included pronounced IFN-like side effects such as pyrexia, chills, myalgia, headache, nausea and malaise, which confirm the relation with the study drug. Importantly, ANA773 was well-tolerated in 13 week toxicology studies at doses that produced robust immune induction (19).

Here, we report the safety profile, pharmacokinetics, antiviral activity and the immunological effects of ANA773 administered orally at four dose levels (800, 1200, 1600 and 2000 mg/day every-other-day) to 34 treatment-naive and treatment-experienced HCV infected patients.

Materials and Methods

Study Design and Organisation

This randomized, placebo-controlled, double-blind, multiple dose escalation, phase 1b

Antiviral activity of ANA773 in chronic HCV patients Final version 22JAN2011

Association Declaration of Helsinki after approval by the institutional review board at each center.

All patients provided written informed consent before participating in any study-related activity.

Initially, HCV patients of any genotype were randomized into three sequential cohorts (800, 1200 and 1600 mg). In the 800, 1200 and 1600 mg groups, 6 patients received oral ANA773 and 2 received placebo every-other-day for 28 days (14 doses). As doses in the first 3 patient dose groups were well tolerated and encouraging antiviral responses were observed at the highest dose level of 1600 mg, a fourth dose group was added to investigate a higher ANA773 dose of 2000 mg. In this dose group, 10 patients (8 active, 2 placebo) were administered 2000 mg of ANA773 or placebo every-other-day for 10 days (5 doses). The shorter treatment period was considered sufficient since immunological and antiviral effects of ANA773 were seen shortly after start of treatment of patients in the previous lower dose groups.

During the first 3 doses of ANA773 (800, 1200 and 1600 mg group), and during all 5 doses in the 2000 mg group, the study was conducted as an inpatient study for careful assessment of safety and tolerability. Study medication (100 mg capsules) and placebo capsules were supplied by Anadys Pharmaceuticals, Inc., San Diego, USA. Patients were allowed to start standard of care (peginterferon + ribavirin) immediately after the study period at the discretion of the patient.

Patients

Key inclusion criteria included male and female patients between 18 to 65 years, with body mass indexes of 18 to 35 kg/m 2 , HCV RNA level ≥ 75 x 10 receiving antiviral therapy or immunomodulatory therapy within 90 days prior to administration of the first dose of study medication were excluded. Patients with chronic stable haemophilia or on stable methadone substitution treatment were eligible.

Safety Assessment

Patients who received at least 1 dose of ANA773 were considered evaluable for safety.

Safety was evaluated by adverse events registration and clinically significant changes from pretreatment baseline in laboratory values, vital signs, electrocardiogram tracings, and findings that were recorded during physical examinations. Dose escalation to the next dose group was based upon the safety and tolerability of ANA773 in all patients in the previous dose group following completion of the 28-day treatment period as determined by the Ethics Committee.

Pharmacokinetic Assessment

Plasma pharmacokinetic samples for assay of the active metabolite of ANA773 and 2 other metabolites were collected at the following times relative to ANA773 dosing: at pre-dose (0 hour), 0.5, 1, 1.5, 2, 3, 4, 8, 12 and 24 hours after dosing on Day 1 in all dose groups. In the 2000 mg group, plasma pharmacokinetic samples were also drawn on Day 9 on the same time points.

Pharmacokinetic analyses were performed by a central laboratory.

Viral Assessment

Siemens INNO-LiPA HCV II assay was used to assess geno(sub)typing of all patients.

Serum samples for HCV RNA analysis were obtained on Day 1 before the first morning dose, followed by sample collection on Day

Statistical Analysis

The determination of the sample size was based on empirical considerations rather than statistical justification. No formal sample size calculations for this study were performed. The sample size of 24 patients in 3 dose groups with 8 patients per dose group (800, 1200 and 1600 mg groups) and an additional 10 patients in the 2000 mg group were considered appropriate for this type of study. Safety data were tabulated and summarized by descriptive statistics.

Results

Demographic and Baseline Characteristics

A total of 54 patients were screened of whom 33 patients were considered eligible. One patient participated twice in the study under two different subject numbers. This was allowed per protocol as the patient received placebo during his first treatment period. The different subject numbers were treated as if they were two individuals and therefore the total number of approved and included subjects was 34. There were no drop outs or withdrawals; all patients completed the study as per protocol. Eighteen patients (53%) were treatment-naive; the remaining 16 patients (47%) had relapsed from a prior full course of IFN-α based therapy. Patients were infected with genotype 1 (n=17), 2 (n=2), 3 (n=10) or 4 (n=5). Of the 34 patients, 30 (88%) were male and 4 (12%) were female. Thirty-one (91%) patients were White/Caucasian. The mean age was 49.

Pharmacokinetics

Pharmacokinetic parameters of the active metabolite of ANA773 are shown in Table 3.

ANA773 was rapidly absorbed in all patients. At all dose levels the active metabolite of ANA773 4 and in Figure 3. The time of nadir for these patients varied between 3 to 9 days after start of treatment and was always 24 h after ANA773 dosing. After reaching nadir, HCV viral load started to increase in several of these patients (Figure 3). Serum viral load levels generally reverted to baseline after cessation of treatment. Viral load reductions of >1.0 log 10 IU/mL during treatment were observed in all genotypes. The pre-treatment viral load ranged from 5.11 to 7.04 log 10 IU/mL in the 7 patients who showed a maximum virological response of >1 log10 IU/mL. A response to ANA773 treatment was seen in both treatment-naive and relapse patients.

The 2 patients with a clinically relevant virological response in the 1600 mg dose group had relapsed from previous treatment, whereas the 5 patients with a clinically relevant virological response in the 2000 mg dose group were all treatment-naive.

Immunological parameters

In total 10 patients had a measurable IFN-α response, defined as having one or more IFNα concentrations above LLQ that increased after ANA773 dosing. For the patients treated with placebo no IFN-α concentrations above LLQ were observed. The incidence of the IFN-α responses appeared to increase with increasing dose of ANA773. The response was rapid: IFN-α concentrations >LLQ were observed at 3 h post-dose in several patients. In most cases, the peak IFN-α concentrations were observed at 6 to 9 h post-dose. The IFN-α concentrations rapidly decreased thereafter and were below LLQ prior to the administration of the next ANA773 dose.

Due to the small number of patients that showed a quantifiable response, it was not clear whether the magnitude of the response increased with dose. Of the 7 subjects who showed a decline in viral load of >1.0 log 10 IU/mL, 4 patients also had a measurable IFN-α response (one patient in the 1600 mg group and 3 patients in the 2000 mg group). However, for 3 of the subjects with a virological response of >1.0 log 10 IU/mL no IFN-α concentrations above LLQ were observed (one patient) or at only one time point (2 patients). In one patient (800 mg group) measurable IFN-α concentrations were present throughout the study period. Since these substantial levels of circulating IFN-α in this subject were present even at baseline, this is considered to be unrelated to the study treatment. Nonetheless, the IFN-α concentrations increased following ANA773 dosing on 4 of the 5 days with post-dose IFN-α assessments; this subject is therefore considered as having a measurable IFN-α response.

The expression of two different IFN-response genes, OAS1 and ISG15, was evaluated as markers of IFN-response. The extent of the increase (mean fold change in mRNA levels) was higher for ISG15 than for OAS1, but other than that the profiles were very similar with regard to timing and dose-response relationship. As previously reported in healthy volunteers dosed with ANA773, the maximum increases in mRNA levels were observed approximately 6 to 12 h post-dose and mean mRNA levels were near baseline by 48 h post-dose. The extent of the induction of these genes appeared to be approximately constant for the first 3-5 dosing occasions.

All 7 patients who showed a decline in viral load of >1.0 log 10 IU/mL, also had clearly increased expression of OAS1 and ISG15.

Mean neopterin concentrations increased for all ANA773 dose levels tested (Figure 4A).

There was a dose-response relationship up to the 1600 mg ANA773 dose. Following the 2000 mg ANA773 dose, the maximum increase in mean neopterin concentrations was slightly lower than following the 1600 mg dose. The highest increases in neopterin concentrations were observed at the beginning of the treatment period. After the first 3 to 4 doses, the neopterin concentrations wained over time such that in most patients neopterin concentrations had returned to values near baseline by the end of the treatment period. In the patients with a viral decline >1.0 log 10 IU/mL a clear increase in neopterin concentrations was observed. At the 2000 mg ANA773 dose level, the 5 patients with a virological response >1.0 log 10 IU/mL had a neopterin mean maximum change (R max, FC ) ranging from 2.02-to 6.48-fold. These patients also had the highest maximum absolute neopterin concentrations. The 3 patients in this dose group with no virological response had a neopterin R max, FC similar to values measured for placebo, indicating no clear increase in neopterin.

Antiviral activity of ANA773 in chronic HCV patients Final version 22JAN2011

An increased 2,5-OAS activity was observed for the 1200 and 1600 mg dose level (Figure 4B). For the 2000 mg dose group no 2,5-OAS activity was determined. Although the mean plots do not indicate a clear dose-response relationship, the mean values for R max, FC did seem to indicate the presence of a dose-response relationship. Since the samples of the 2000 mg dose group were not analyzed, only a limited evaluation of a possible relation between increased 2,5-OAS activity and virological response could be done. The 2 patients who had a decrease in viral load of >1.0 log 10 IU/mL, did not appear to have the highest increases or absolute values of 2,5-OAS activity in their dose groups. The increase in mean 2,5-OAS activity was observed starting from 24 h after the first dose. A further increase was observed after the second dose;

thereafter the 2,5-OAS activity appeared to remain constant throughout the treatment period.

Discussion

The present study is the first clinical trial to evaluate ANA773, an oral inducer of endogenous interferons that acts via TLR7, in chronic hepatitis C patients. This multiple dose escalation study was designed to explore the safety, tolerability, pharmacokinetics, pharmacodynamics and antiviral activity of ANA773.

Our primary objective was to investigate safety and tolerability of ANA773. There were no serious AEs observed and there were no premature discontinuations. The most frequently reported AEs, which were considered to be related to the study medication, were flu-like symptoms, which were only observed early following the start of treatment. Since these symptoms are known IFN-like side effects, the observed AE profile is consistent with the expected mechanism of action of the compound. With increasing dose of ANA773, a trend was seen towards a higher frequency and stronger intensity of flu-like symptoms. Overall, treatment with multiple doses of ANA773 at dose levels up to 2000 mg were well tolerated by the 26 patients chronically infected with HCV. The same mild profile of side effects was seen with an intravenous TLR7 agonist (isatoribine), although only few flu-like symptoms were seen in this study [START_REF] Ghany | Diagnosis, management, and treatment of hepatitis C: an update[END_REF]. In contrast, oral administration of resiquimod (a TLR7 and TLR8 agonist) was associated with typical interferon-related adverse events which limited dose escalation [START_REF] Horsmans | Isatoribine, an agonist of TLR7, reduces plasma virus concentration in chronic hepatitis C infection[END_REF].

ANA773 was rapidly absorbed following oral administration, and all three metabolites were detected in plasma at the first sampling time point 30 minutes after start of treatment. Systemic exposure to the active metabolite of ANA773 was dose proportional, although the C max increased with dose in a less than dose proportional manner. It is possible that the rate of conversion of ANA773 into the active metabolite may become more limiting at the higher doses. An additional factor contributing to the less than dose proportional enhancement of the C max of the active metabolite of ANA773 could be the high number of capsules administered which may limit the efficiency of the dissolution and absorption process. For the C max of metabolite A no deviations from dose proportionality were observed, which would also be expected if dissolution and absorption was a major factor. The pharmacokinetic profile of ANA773 in chronic HCV infected patients was similar to healthy volunteers (data not shown). For the 2000 mg group it was shown that all three metabolites were below LLQ at pre-dose on Day 9, indicating that no accumulation during multiple dosing occurred.

Repeated administration with 1600 mg or 2000 mg of ANA773 resulted in a mean maximum viral load decrease of -0.97 and -1.26 log 10 IU/mL, respectively, compared to -0.34 log 10 IU/mL in placebo-treated HCV patients (p=0.037 at 2000 mg relative to the placebo group). At the 1600 mg dose level 2 of 6 HCV patients and at the 2000 mg dose level 5 of 8 HCV patients had a decline in viral load of >1.0 log 10 IU/mL. The decline of serum HCV RNA levels was generally observed within the first several days, and HCV RNA levels remained reduced until the end of therapy in most patients. Interestingly, HCV viral load decline did not resemble the typical biphasic pattern as seen with interferon-based regimens, which may be due to a different antiviral mechanism, suboptimal drug dosing, or due to the specific dosing regimen.

Antiviral activity of ANA773 in chronic HCV patients Final version 22JAN2011

Stimulation with TLR7 agonists induces the production of type I interferons, proinflammatory cytokines and the expression of co-stimulatory molecules on various cell populations [START_REF] Poordad | Boceprevir for untreated chronic HCV genotype 1 infection[END_REF][START_REF] Hemmi | Small antiviral compounds activate immune cells via the TLR7 MyD88-dependent signaling pathway[END_REF][START_REF] Savage | A phase I clinical trial of imiquimod, an oral interferon inducer, administered daily[END_REF]. In all patients with a decline of viral load of more than 1.0 log 10 these factors were induced or activated, demonstrating the immunomodulatory activity of ANA773 leading to antiviral effects. Our findings demonstrated that the decline in serum HCV RNA levels correlated with exposure to the active metabolite of ANA773 as well as neopterin levels, respectively (Figure 5).

Further studies are needed to examine why some patients respond to ANA773 as demonstrated by activation of innate immune responses, whereas no decline of viral load is observed in these patients. We observed that HCV genotype and the baseline HCV RNA levels were not predictive of the responsiveness to administration of the TLR7 agonist.

Whereas current peginterferon administration consists of a single IFN subtype (2a or 2b), TLR activation induces several IFN subtypes which may potentially improve antiviral efficacy, but may also lead to more cytokine mediated adverse events; however, both improved antiviral efficacy as deterioration of the AE profile were not seen (with respect to peginterferon), suggesting that the optimal dosage of ANA773 remains to be determined.

In this proof of concept study conclusions are limited by the number of subjects and the heterogeneity of the study population. But due to its heterogenous population it is possible to conclude that immune activation and viral decline can be achieved by ANA773 in HCV infected patients, irrespective of genotype, viral load, sex or previous response to IFN-based therapies.

In conclusion, oral administration of ANA773 was safe and modest anti-HCV activity was seen in combination with marked therapy induced activation of the immune system. Potential advantages of ANA773 relative to weekly subcutaneous injections with peginterferon include the route of administration (oral versus injection) and a good safety profile. For the near future it is likely that peginterferon will remain the cornerstone of antiviral combination therapies. However, with the advent of direct antiviral compounds targeting the HCV protease and polymerase proteins and leading to strong viral suppression, the clinical development of immune activators, such as TLR-agonists, may be advantageous. In this study, we show that forced viral decline combined with activation of the immune system by the TLR7 agonist ANA773 is safe and exhibits modest anti-HCV activity. In addition, our findings pave the way for future studies to explore whether higher doses of ANA773 administered via the oral route will further enhance the antiviral efficacy of this TLR7 agonist in lowering the serum HCV RNA levels in chronic HCV patients. It will be of particular interest to investigate this drug in combination with direct antivirals or with ribavirin, as currently is being planned. (24 hours after the last dose). Patients in the 1600 mg group had samples taken pre-treatment, on days 1-6, and then on days 13, 21, 27 and 28 (24 hours after the last dose). For these patients, the period between days 6-28 is represented by a dashed line. Patient genotypes (GT): green; GT3, orange; GT1A, light blue; GT2, black; GT3, dark blue; GT4, red; GT4, purple; GT3.

F

  Pharmacokinetic parameters were summarized by descriptive statistics at each dose level and by study day. Analysis of variance (ANOVA) was used to compare the dose normalized C max and AUC values across dose groups. The log transformed HCV RNA values were summarized by descriptive statistics according to dose group and study. ANOVA was used to compare the difference between each time point and baseline, the nadir and baseline, and between the end-oftreatment and baseline across dose groups (including placebo). Immunological parameters for serum cytokines (IFN-α), RNA (OAS and ISG15), neopterin and 2,5-OAS activity wereAntiviral activity of ANA773 in chronic HCV patients Final version 22JAN2011summarized by descriptive statistics according to dose group and study day. The software used for all summary statistics and statistical analyses were SAS version 9.1.3 (SAS Institute, Inc.) and comparable desktop parameters. For the calculation of some pharmacokinetic parameters WinNonlin (Pharsight, Inc.) version 5.0.1 was used. For the calculation of correlation between variables the Spearman non-parametric test was used. Any p-values calculated as a result of the statistical analysis were interpreted in accordance with the exploratory nature of this study.

Figure 2 .

 2 Figure 2. Individual viral load data of the 7 patients with >1.0 log 10 IU/mL decline during treatment are shown in Table4and in Figure3. The time of nadir for these patients varied between 3 to 9

Figure 1 .

 1 Figure 1. Mean plasma concentration of the active metabolite of ANA773 on Day 1 (all groups) and Day 9 (2000 mg only) SD: single dose; MD: multiple dose.

Figure 2 .

 2 Figure 2. Mean maximum HCV RNA response of all patients during treatment with ANA773 Mean (+ SEM) viral load nadir values are plotted for each dose group. The maximum viral load decline at 2000 mg is statistically significant relative to the placebo group. The number of patients in each group with a maximum decline in viral load of >1.0 log10 IU/mL is displayed above the plot. NOTE: patients in the 800, 1200 and 1600 mg group were dosed for 28 days. Patients in the 2000 mg group were dosed for 10 days.

Figure 3 .

 3 Figure 3. Viral load time course in patients with >1.0log 10 decline during treatment Five 2000 mg patients (closed circles) and two 1600 mg patients (open circles) had a maximum decline in viral load of >1 log10 during treatment.Each line represents the viral load change for one of these patients during the treatment period (10 days for the 2000 mg group, 28 days for the 1600 mg group). Patients in the 2000 mg group had samples taken pre-treatment (PT) and on days 1 through 10

Figure 4A .

 4A Figure 4A. Mean neopterin concentration throughout the study NOTE: patients in the 800, 1200 and 1600 mg group were dosed for 28 days. Patients in the 2000 mg group were dosed for 10 days.

Figure 4B .Figure 5 .

 4B5 Figure 4B. Mean 2,5-OAS activity throughout the study

  study was conducted at three sites in The Netherlands from October 2008 until August 2009. This study was conducted in accordance with Good Clinical Practice and with the World Medical

	F o r
	P
	e
	e r
	R
	e
	v i e
	w

Antiviral activity of ANA773 in chronic HCV patients Final version 22JAN2011

	F o r
	P
	e
	e r
	R
	3 IU/mL, and clinical laboratory e
	evaluations consistent with chronic hepatitis C infection as defined by the protocol. Treatment-naive or relapse patients were allowed. Relapse was defined as undetectable HCV RNA at completion of a previous IFN-based treatment, but positive HCV RNA during follow-up. Key v i e
	exclusion criteria included decompensated liver disease, findings consistent with Child Pugh B/C liver cirrhosis, and co-infection with HIV or HBV. Previous non-responders to IFN-based therapies w
	Comment [W5]: see JB5
	were excluded due to expected unresponsiveness to immunomodulatory therapy. Patients

Antiviral activity of ANA773 in chronic HCV patients Final version 22JAN2011 Immunological parameters Cytokine

  serum samples (IFN-α) were obtained on Day 1, 2, 3, 5, 13, 27, 28, and 34 (800, 1200 and 1600 mg groups). In the 2000 mg group, cytokine serum samples were obtained on Day 1, 3, 5, 7 and 9. Samples were analyzed by Enzyme-Linked ImmunoSorbent Assay (ELISA) at Alta Analytical Laboratory, San Diego, USA. Serum RNA samples (OAS1 and ISG15) were obtained on Day 1, 2, 3, 4, 5, 6, 13, 21, 27, 28, 34 (800, 1200 and 1600 mg groups) and on Day 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 16 (2000 mg group) by Anadys Pharmaceuticals, San Diego, USA, using branched DNA (bDNA) assay. Neopterin and 2,5-OAS activity samples were collected at Day 1, 2, 3, 4, 5, 6, 13, 21, 27, 28, 34, 41 (800, 1200 and 1600 mg groups) and at Day 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 16 (2000 mg group). Neopterin samples were analyzed by ELISA at Alta Analytical

	F o r
	P
	e
	e r
	R
	e
	2, 3, 4, 5, 6, 13, 21, 27, 28, 34 and 41 (800, 1200 and 1600 mg groups) or on Day 2, 3, 4, 5, 6, 7, 8, 9, 10 and 16 (2000 mg group). HCV RNA measurements v i e
	were performed by a central laboratory using the Roche COBAS AmpliPrep/COBAS TaqMan HCV Test with a lower limit of quantification (LLQ) of 43 IU/ml. A >1.0 log 10 reduction in viral load during w
	treatment was considered a clinically relevant decline.

Laboratory, San Diego, USA. 2,5-OAS activity was analyzed by radioimmunoassay (RIA) at PRA International, Assen, The Netherlands.

Antiviral activity of ANA773 in chronic HCV patients Final version 22JAN2011 Safety Safety

  data including serious adverse events, premature discontinuations and hematologic effects are shown in Table2. Overall, multiple oral doses of ANA773 were safe and generally well tolerated at all dose levels tested up to 2000 mg. There were no deaths or other SAEs during the study. The majority of AEs were of mild or moderate intensity. The total number of AEs was highest in the 1600 mg (66 AEs) and 2000 mg (58 AEs) dose group. The number of AEs reported by the placebo-treated patients was only slightly less (49 AEs). A total of 167 AEs were considered possibly or probably related to the study medication (71% of the total number of AEs).The most frequently reported related AEs were headache, myalgia, fatigue, abdominal pain and dizziness. The occurrence of headache appeared to increase with dose, whereas for the other frequently reported related AEs a relation with dose was not observed. In addition, nausea (total of

	F o r
	P
	e
	e r
	R
	e
	v i e
	w

5 years (range, 32-64 years) and mean BMI was 26.35 kg/m 2 (range, 18.3-34.6 kg/m 2 ).

Demographic and other baseline characteristics of the randomized patients are shown in Table

1

.

6 events) considered related to the study medication occurred in the highest dose groups only. At the 1600 mg dose level 2 patients experienced flu-like symptoms (e.g. headache, myalgia, chills, malaise) at the beginning of their dosing period. However, these symptoms rapidly disappeared during continued dosing, even when they were initially of moderate to severe intensity. At the 2000 mg dose level 3 patients treated with ANA773 experienced flu-like symptoms, mainly at the first 2 dosing occasions. One of these patients started to use paracetamol as prophylactic treatment and the remaining doses were well tolerated. Individual dose reductions, as were allowed per protocol, were not necessary. Trends for decreases in leukocytes (nadir -4.06 x 10 9 /l on Day 6), thrombocytes (nadir -38.4 x 10 9 /l on Day 6) and lymphocytes (nadir -1.23 x 10 9 /l on Day 2) were observed for the highest dose level (2000 mg ANA773), but were not considered clinically significant. No other clinically significant abnormalities in vital signs, ECG, clinical laboratory values and physical examination were observed.

Antiviral activity of ANA773 in chronic HCV patientsTable 1 .

 1 Demographics

					21
	800 mg	1200 mg	1600 mg	2000 mg	Placebo
	(n=6)	(n=6)	(n=6)	(n=8)	(n=8)
	o r F o r F o r F o r Male, n (%) F o r
		P P P P P
		e e e e e
		e r e r e r e r e r
			R R R R R
			e e e e e
			v i e v i e v i e v i e v i e
				w w w w w
	Antiviral activity of ANA773 in chronic HCV patients			Final version 22JAN2011

Antiviral activity of ANA773 in chronic HCV patientsTable 2 .

 2 Safety results

			Total	800 mg	1200 mg	1600 mg	2000 mg	Placebo
			(n=34)	(n=6)	(n=6)	(n=6)	(n=8)	(n=8)
	Serious adverse events (n)	0
	Discontinuations (n)	0
	Adverse events (n)	235	31	31	66	58	49
	Patients with AE (%) 1
	-	Headache	53	33	33	67	75	55
	-	Fatigue	41	50	33	67	38	25
	-	Myalgia	38	33	17	50	50	38
	-	Abdominal pain	32	17	17	50	13	63
	-	Dizziness	24	50	17	50	0	13
	-	Back pain	21	0	17	50	13	25
	----	Constipation Nausea Dyspepsia Malaise	18 18 15 15 F o 0 0 0 17 r	0 0 0 0	50 33 33 33	0 50 25 25	38 0 13 0
	-	Chills	12	0	17	33	13	0
	---	Diarrhoea Pain in extremity Pyrexia	12 12 12	17 17 P 0 e 0 0 0	17 17 0	0 13 38	25 13 13
	Hematologic abnormalities -Anemia 2 -Leukocytopenia 3	0 0		e r
	-1 AEs listed that occurred in >10% of all patients Trombocytopenia 4 0			R
	2 Hemoglobin level < 6.2 mmol/l 3 Leukocyte count < 1.5x10 9 /l 4 Platelet count < 50x10 9 /l				e v i e
							w

Antiviral activity of ANA773 in chronic HCV patients

Table 3 .

 3 Pharmacokinetic parameters of plasma ANA773 active metabolite

	Dose (mg)	Day	Cmax (ug/ml)	Tmax (h)	AUCinf (ug.h/ml)	T1/2 (h)
	800	1	10.69 (7.57-15.3)	0.509 (0.50-1.00)	24.96 (21.9-29.1)	3.254 (2.09-3.96)
	1200	1	11.99 (8.10-16.9)	1.000 (0.50-1.50)	31.53 (26.2-34.7)	2.901 (1.87-3.75)
	1600	1	13.94 (9.17-19.0)	1.000 (0.50-2.00)	40.76 (33.3-49.4)	3.375 (2.66-3.85)
	2000	1	15.87 (10.9-23.9)	1.000 (0.50-1.50)	51.06 (36.7-69.4)	2.995 (1.55-4.29)
	2000	9	15.50 (13.1-19.3)	1.000 (0.50-1.50)	46.42 (35.0-60.4)	3.237 (2.67-3.65)
	For all parameters except tmax the geometric mean (range) is presented; for tmax the median (range) is presented.

Antiviral activity of ANA773 in chronic HCV patientsTable 4 .

 4 Individual HCV RNA data for patients with a maximal viral decline of >1.0 log 10 IU/mL Mean of Day -1 and pre-dose Day 1; 2 Maximum negative change from baseline value at Day 28 or at Day 10; 3 Average of the Day 27 and Day 28 or Day 9 and Day 10 (2000 mg only) HCV RNA concentrations.

	Subject	HCV	Sex	Baseline	Nadir	Day (N)	Day 1	Day 5	Day 28	EOT /	Final
		genotype		viral			(N)	(N)	(N) / Day	EOT9/10 3	follow-up
				load 1					10 (N) 2
	1600 mg
	Patient 1	3	M	7.04	-1.15	4	-0.74	-0.70	-0.49	6.52	6.86
	Patient 2	1A	F	6.59	-2.52	6	-1.64	-2.52	-2.24	4.41	6.02
	2000 mg
	Patient 3	3A	F	6.77	-1.63	10	-1.62	-1.30	-1.63	5.35	6.10
	Patient 4	4D	M	6.23	-1.35	6	-1.01	-1.35	-0.55	6.10	5.56
	Patient 5	4A	M	5.11	-1.53	4	-0.65	-1.12	-0.28	4.72	4.67
	Patient 6	3A	M	6.99	-3.10	4	-2.16	-3.07	-1.68	5.25	5.91
	Patient 7	2	M	6.75	-1.78	8	-0.55	-1.03	-1.07	5.75	6.95
				F o r
						P
						e
						e r
								R
								e
								v i e
										w

1

Antiviral activity of ANA773 in chronic HCV patients

(100) 6 (100) 4 (67)

(88)

Acknowledgement

The authors thank the patients who generously agreed to participate in this clinical study.

We also acknowledge the contribution of the following individuals: C. van der Ent, I. Brings (Clinical Research Bureau, Erasmus MC, Rotterdam), M.W. Peters (study nurse at Academic Medical Center, Amsterdam), B. Eam, M.V. Sergeeva and T.W. Harding (Anadys Pharmaceuticals, San Diego).

This study was sponsored by Anadys Pharmaceuticals, Inc.

Squibb and receives grant/research support from Roche, Novartis, Gilead, Schering-Plough, and Bristol-Myers Squibb. R.J. de Knegt serves as advisor for Schering-Plough, receives grant/research support from Schering-Plough, Roche, and Medtronic, and speaks and teaches for Schering-Plough. H.W. Reesink serves as advisor for Roche Molecular Diagnostics, Anadys, , Arrows, and Gilead, consults for PRA-International, Tibotec, GlaxoSmithKline, Chiron Novartis, and Roche Therapeutics, and receives grant/research support from Schering-Plough,

Antiviral activity of ANA773 in chronic HCV patients

Final version 22JAN2011


Jan Bulla

Christophe Chesneau
email: ches-neau@math.unicaen.fr

Maher Kachour
email: maher.kachour@univ-nantes.fr

A bivariate first-order signed integer-valued autoregressive process

Keywords: Bivariate integer-valued time series, INAR models, SINAR models, Bivariate Skellam distribution. 2000 Mathematics Subject Classification 62M10, 62M20

Bivariate integer-valued time series occur in many areas, such as finance, epidemiology, business etc. In this paper, we present bivariate autoregressive integer-valued time series models, based on the signed thinning operator. Compared to classical bivariate INAR models, the new processes have the advantage to allow for negative values for the time series and the autocorrelation functions. Strict stationarity and ergodicity of the processes are established. The moments and the autocovariance functions are determined. Some methods for estimating the model parameters are considered and the asymptotic properties of the obtained estimators are derived. Simulation experiments as well as analysis of real data sets are carried out to assess the models' performance.

Introduction

In many practical situations we have to deal with integer-valued times series. Because of the integer character of the observed series, a standard autoregressive process cannot be considered. Many attempts have therefore been made to define univariate models having the same properties as a real-valued AR process and respecting the integer nature of the observations. Some important contributions are due to [START_REF] Mckenzie | Some simple models for discrete variate time series[END_REF], [START_REF] Al-Zaid | First-order integer-valued autoregressive (INAR(1)) process[END_REF] and [START_REF] Du | The integer-valued autoregressive (INAR(p)) model[END_REF]. They used the binomial thinning operator of [START_REF] Steutel | Discrete analogues of self-decomposability and stability[END_REF] to define stochastic processes, denoted by INAR, which resemble AR(1) or AR(p) processes. Gauthier and Latour (1991) also considered a more general version called the GINAR(p) model, based on the generalized Steutel and van Harn operator. For a recent review on discrete-valued time series models, we refer to, e.g., [START_REF] Kedem | Regression models for time series analysis[END_REF], [START_REF] Mckenzie | Discrete variate time series[END_REF], [START_REF] Jung | Binomial thinning models for integer time series[END_REF] and [START_REF] Weiß | Thinning operations for modeling time series of counts a survey[END_REF].

The literature on multivariate integer-valued time series models is limited. A first model (by [START_REF] Franke | Multivariate first-order integer-valued autoregressions[END_REF], which introduced the multivariate INAR(1) process, is based on the classical binomial thinning operator. [START_REF] Latour | The multivariate GINAR(p) process[END_REF] introduced a multivariate GINAR(p) model, based on the generalized thinning operator. Recently, there have been more attempts to derive bivariate integer-valued autoregressive type models. The interested reader is referred to Brannas and Nordstrom (2000), [START_REF] Quoreshi | Bivariate Time Series Modeling of Financial Count Data[END_REF], [START_REF] Heinen | Multivariate autoregressive modeling of time series count data using copulas[END_REF] and [START_REF] Pedeli | A bivariate Poisson INAR(1) model with application[END_REF].

All models mentioned above (univariate and multivariate) are subject to the same intrinsic limitations. By construction, they cannot fit a time series with negative observations. Moreover, since regression coefficients must be positive, the modelling of series with possible negative correlations is not possible. In order to avoid the shortcomings of previous models, [START_REF] Kachour | A p-order signed integer-valued autoregressive (SINAR(1)) model[END_REF] introduced a general univariate class, denoted by SINAR, based on a modified version of the generalized thinning operator, called the signed thinning operator. Indeed, SINAR allows for negative values both for the series and its autocorrelation function, and can be seen as an extension of the GINAR model of Gauthier and Latour (1991) on Z.

The purpose of this paper is to introduce a simple bivariate integer-valued autoregressive model, based on the signed thinning operator, denoted by B-SINAR(1). More precisely, we present two model specifications: First, with full autoregressive matrix and independent noise components. Second, with diagonal autoregressive matrix and correlated noise. Thus, the first specification is considered as an extension of the multivariate GINAR model on Z 2 (in the bivariate case). Moreover, the second specification can be seen as a general version of the model presented by [START_REF] Pedeli | A bivariate Poisson INAR(1) model with application[END_REF].

The remainder of this paper is organized as follows. In Section 2 we present our main results. More precisely, in Section 2.1, two specifications of the first-order bivariate signed integer-valued autoregressive model are introduced. Moreover, we provide theoretical results about the stationarity and the moments of both B-SINAR(1) processes. In Section 3, the estimation of the model parameters is considered. Simulation experiments as well as analysis of real data are discussed in Sections 4 and 5. The proofs of all results are postponed to Section 6.

Main results

The B-SINAR(1) process

The so-called signed thinning operator, originally proposed by [START_REF] Latour | An integer-valued bilinear type model[END_REF], is a natural extension of the [START_REF] Steutel | Discrete analogues of self-decomposability and stability[END_REF] operator to Z-valued random variables and is defined as follows.

Definition 1 (Signed thinning operator ) Let (Y i ) i∈Z be a sequence of i.i.d. integer-valued random variables with F as common distribution, independent of an integer-valued random variable X. The signed thinning operator, denoted by F •, is defined by

F • X =            sign(X) |X| i=1 Y i , if X = 0, 0 otherwise, (1)
where, for an integer x, sign(x) = 1 if x > 0 and -1 if x < 0. The sequence (Y i ) i∈Z is referred to as a counting sequence.

The model introduced in this paper is based on a random matrix operator defined below.

Definition 2 (Signed matrix thinning operator ) Let (Z 1 , Z 2 ) τ be an integervalued random vector. The signed matrix thinning operator, denoted by F , is defined by

F       Z 1 Z 2       =       F 1,1 • Z 1 + F 1,2 • Z 2 F 2,1 • Z 1 + F 2,2 • Z 2      
, where -F i,j is the common distribution of the i.i.d. counting sequence (Y (i,j) k

) k∈Z , all the counting sequences associated to the operators F i,j • are mutually independent for any (i, j) ∈ {1, 2} 2 .

Definition 3 (B-SINAR(1) process of first kind ) A bivariate process (X t,1 , X t,2 ) t∈Z is said to be a B-SINAR(1) (for Bivariate Signed INteger-valued AutoRegressive) process of first kind if it admits the following representation

      X t,1 X t,2       = F       X t-1,1 X t-1,2       +       ε t,1 ε t,2       , (2)
where (ε t,j ) t∈Z is an independent stationary process such that support(ε 1,j ) ⊆ Z, it is independent of all counting sequences of the model, -(ε t,1 ) t∈Z and (ε t,2 ) t∈Z are independent for any j ∈ {1, 2}.

Adopting the notations of Definitions 2 and 3, we set

γ i,j = E(Y (i,j) 1
),

β i,j = V(Y (i,j) 1
),

µ j = E(ε 1,j ), σ 2 εj = V(ε 1,j )
and

M =       γ 1,1 γ 1,2 γ 2,1 γ 2,2       (3)
for any (i, j) ∈ {1, 2} 2 .

The following theorem provides the conditions which ensure the ergodicity and the stationarity of the B-SINAR(1) process of first kind.

Theorem 1 Consider a B-SINAR(1) model of first kind, as defined in (2). Assume that:

1. ε 1 = (ε 1,1 , ε 1,2 ) charges all points of E = Z 2 , i.e., P(ε 1 = x) > 0, ∀ x ∈ Z 2 , 2. ρ(M )
, the spectral radius of the autoregressive matrix of the process M (3), is smaller than 1.

Then, the process has a unique stationary solution, denoted by X = (X 1 , X 2 ), and it is also ergodic. Moreover, E(X 2 j ) < ∞ for any j ∈ {1, 2}.

In the following, we derive some properties of the B-SINAR(1) process of first kind under the assumptions of Theorem 1. Let

X 1 d = X t,1 and X 2 d = X t,2
. By the stochastic representation in (2), we have

E(X 1 ) = µ 1 (1 -γ 2,2 ) + γ 1,2 µ 2 (1 -γ 1,1 )(1 -γ 2,2 ) -γ 1,2 γ 2,1 , V(X 1 ) = β 1,1 E(|X 1 |) + β 1,2 E(|X 2 |) + γ 2 1,2 V(X 2 ) + σ 2 ε1 1 -γ 2 1,1 , E(X 2 ) = µ 2 (1 -γ 1,1 ) + γ 2,1 µ 1 (1 -γ 1,1 )(1 -γ 2,2 ) -γ 1,2 γ 2,1 and V(X 2 ) = β 2,1 E(|X 1 |) + β 2,2 E(|X 2 |) + γ 2 2,1 V(X 1 ) + σ 2 ε2 1 -γ 2 2,2
.

Moreover, the autocorrelation functions are given by

ρ X1 (k) = γ k 1,1 + γ 1,2 V(X 1 ) k j=1 γ k-j 1,1 Cov (X t,1 , X t+j-1,2 ) and ρ X2 (k) = γ k 2,2 + γ 2,1 V(X 2 ) k j=1 γ k-j 2,2 Cov (X t,2 , X t+j-1,1 ) ,
for any k ∈ N * , where

Cov (X t,1 , X t+h,2 ) = γ 2,1 h i=1 γ h-i 2,2 Cov (X t,1 , X t+i-1,1 ) + γ h 2,2 1 -γ 1,1 γ 2,2 -γ 1,2 γ 2,1 (γ 1,1 γ 2,1 V(X 1 ) + γ 1,2 γ 2,2 V(X 2 )) , Cov (X t,2 , X t+h,1 ) = γ 1,2 h i=1 γ h-i 1,1 Cov (X t,2 , X t+i-1,2 ) + γ h 1,1 1 -γ 1,1 γ 2,2 -γ 1,2 γ 2,1 (γ 1,1 γ 2,1 V(X 1 ) + γ 1,2 γ 2,2 V(X 2 ))
for any h ∈ N.

Remark 1 Compared to all bivariate integer-valued autoregressive models based on classical or generalized thinning operators, the B-SINAR(1) process of first kind has several advantages. More precisely, our process is constructed in a way such that it can be used to analyze time series with negative values. Moreover, the second assumption in Theorem 1 is similar to the corresponding condition ensuring stationarity of a real-valued bivariate AR(1) process. Therefore, the B-SINAR(1) process of first kind can reproduce auto-and cross-correlation functions with negative values.

Particular case of the B-SINAR(1) process of first kind

Next, we present a modified version of the B-SINAR(1) process of first kind.

Definition 4 (B-SINAR(1) process of second kind ) A bivariate process (X t,1 , X t,2 ) t∈Z is said to be a B-SINAR(1) process of second kind if each component admits the same representation of a SINAR(1) process, i.e.,

X t,1 = F 1 • X t-1,1 + ε t,1 , X t,2 = F 2 • X t-1,2 + ε t,2 , (4)
where

-F j is the common distribution of the i.i.d. counting sequence (Y (j)
k ) k∈Z for any j ∈ {1, 2}, -(ε t,j ) t∈Z is an independent stationary process such that support(ε 1,j ) ⊆ Z, it is independent of all counting sequences of the model for any j ∈ {1, 2}.

Adopting the notations in Definition 4, we set

γ j = E(Y (j) 1 ), β j = V(Y (j) 1 )
for any j ∈ {1, 2}

Remark 2 In Definition 4, it is understood that ε t,1 and ε t,2 are not independent, which implies Cov(ε 1,1 , ε 1,2 ) = 0.

Remark 3 The model defined in (4) can be seen as a particular case of the above B-SINAR(1) process of first kind, where the signed matrix thinning operator equals

F =       F 1 • 0 0 F 2 •      
and thus the autoregressive matrix has the following form

M =       γ 1 0 0 γ 2       . (5)
Hence, using Theorem 1, the following corollary gives the conditions which ensure the ergodicity and the stationarity of the B-SINAR(1) process of second kind.

Corollary 1 Consider a B-SINAR(1) model of second kind, as defined in (4).

Assume that:

1. ε 1 = (ε 1,1 , ε 1,2 ) charges all points of E = Z 2 , i.e., P(ε t = x) > 0, ∀ x ∈ Z 2 , 2. max(|γ 1 |, |γ 2 |) < 1.
Then, the process has a unique stationary solution, denoted by X = (X 1 , X 2 ), and it is also ergodic. Moreover, E(X 2 j ) < ∞ for any j ∈ {1, 2}. Hence, under second order stationarity conditions we obtain

E(X j ) = µ j 1 -γ j , (6)
V(X j ) =

β j E(|X j |) + σ 2 εj 1 -γ 2 j (7) and ρ Xj (k) = γ k j , ∀ k ∈ N (8)
for any j ∈ {1, 2}. Furthermore, we have

Cov (X t,1 , X t+h,2 ) = γ h 2 1 -γ 1 γ 2 Cov (ε t,1 , ε t,2 ) (9)
and

Cov (X t,2 , X t+h,1 ) = γ h 1 1 -γ 1 γ 2 Cov (ε t,1 , ε t,2 ) (10)
for any h ∈ N.

Remark 4 Since the autoregressive matrix defined in ( 5) is diagonal, one can see that correlation between the innovations is the only source of dependence between the two series (X t,1 ) t∈Z and (X t,2 ) t∈Z .

Remark 5 This model can be also viewed as a Z 2 -extension of the one studied in [START_REF] Pedeli | A bivariate Poisson INAR(1) model with application[END_REF]. In addition, due to the second assumption of Corollary 1, the regression coefficients of our process can have negative values, which leads to possibly negative values of the auto-and cross-correlation functions (which cannot be achieved by the model of [START_REF] Pedeli | A bivariate Poisson INAR(1) model with application[END_REF] having coefficients belonging to [0, 1]).

3 Parameter estimation

For B-SINAR(1) process of first kind

Next, we state asymptotic results for the estimation of the parameters of B-SINAR(1) process of first kind, defined by (2). In more detail, we consider the Conditional Least Squares (CLS) estimator to determine γ i,j and µ j , i.e., the means of distributions F i,j and ε j , respectively, for any (i, j) ∈ {1, 2} 2 . Let

θ =       γ 1,1 γ 1,2 µ 1 γ 2,1 γ 2,2 µ 2       .
We suppose that θ belongs to a compact parameter space Θ. Let X 0 , X 1 , • • • , X n be a bivariate sample from B-SINAR(1) process of first kind, where X i = (X i,1 , X i,2 ) for any i ∈ {0, 1, . . . , n}. Thus, the CLS estimator equals

θCLS =       γ1,1 γ1,2 μ1 γ2,1 γ2,2 μ2       = argmin θ∈Θ S n (θ), with S n (θ) = n t=1 X t -θY t-1 2 ,
where • denotes the Euclidian norm on R 2 and Y t-1 = (X t-1,1 , X t-1,2 , 1) τ . Note that we have the following explicit expressions θCLS

1 =            γ1,1 γ1,2 μ1            = Q -1 q 1 , and θCLS 2 =            γ2,1 γ2,2 μ2            = Q -1 q 2 , (11)
where

Q =            n t=1 X 2 t-1,1 n t=1 X t-1,1 X t-1,2 n t=1 X t-1,1 n t=1 X t-1,1 X t-1,2 n t=1 X 2 t-1,2 n t=1 X t-1,2 n t=1 X t-1,1 n t=1 X t-1,2 n            , q 1 =            n t=1 X t,1 X t-1,1 n t=1 X t,1 X t-1,2 n t=1 X t,1            , and q 2 =            n t=1 X t,2 X t-1,1 n t=1 X t,2 X t-1,2 n t=1 X t,2           
Theorem 2 Consider a B-SINAR(1) process of first kind as defined in (2). Suppose that the assumptions of Theorem 1 are satisfied and that E(|X j | 4 ) < ∞ for any j ∈ {1, 2}. Then, we obtain √ n θCLS

1 -θ * 1 L -→ N 0, V -1 ΣV -1 and √ n θCLS 2 -θ * 2 L -→ N 0, V -1 Γ V -1 , where V = lim n→∞ 1 n Q = E Y t-1 Y τ t-1 , Σ = E Y t-1 Y τ t-1 (X t,1 -γ 1,1 X t-1,1 -γ 1,2 X t-1,2 -µ 1 ) 2 and Γ = E Y t-1 Y τ t-1 (X t,2 -γ 2,1 X t-1,1 -γ 2,2 X t-1,2 -µ 2 ) 2
for the CLS estimators defined by (11).

Remark 6 According to Theorem 1 of [START_REF] Kachour | A p-order signed integer-valued autoregressive (SINAR(1)) model[END_REF], one can deduce that E(|X j | 4 ) < ∞ if F i,j and ε 1,j have a fourth moment for any (i, j) ∈ {1, 2} 2 .

For B-SINAR(1) process of second kind

In this section we discuss the estimation problem of a B-SINAR(1) process of second kind under some additional parametric assumptions. Indeed, we assume that F j , the common distribution of the i.i.d. counting sequence (Y

(j) k ) k∈Z , is given by P(Y (j) 1 = -1) = (1-α j ) 2 , P(Y (j) 1 = 0) = 2α j (1-α j ), P(Y (j) 1 = 1) = α 2 j , (12)
with α j ∈ (0, 1) for any j ∈ {1, 2}. In other words, we set Y (j) 1 d = Z j -1 for any j ∈ {1, 2}, where Z j ∼ Binomial(2, α j ).

Remark 7 From (12), for any j ∈ {1, 2}, we have

|γ j | = |E(Y (j) 1 )| = |2α j -1| < 1, β j = V(Y (j) 1 ) = 2α j (1 -α j ). Remark 8 Let j ∈ {1, 2}, x ∈ Z * and y ∈ Z. Than holds P (F j • x = y) = P T (j) x = sign(x)(x + y) ,
where T (j)

x ∼ Binomial(2|x|, α j ), i.e.,

P (F j • x = y) = 2|x| |x| + y sign(x) α |x|+y sign(x) j (1-α j ) |x|-y sign(x) 1 {y∈{-|x|,...,|x|}} ,
with 1 denoting the indicator function.

Remark 9 In fact, the distribution F j , defined by ( 12), can be seen as a natural extension of the binomial distribution, which is used in the definition of the classical thinning operator.

Moreover, we assume that the innovations of the two series in (4) jointly follow the bivariate Skellam distribution presented below.

Definition 5 (BSkellam(λ 0 , λ 1 , λ 2 ) distribution) Let λ 0 ≥ 0, λ 1 > 0 and λ 2 > 0. We say that the bivariate random variable (ε 1 , ε 2 ) has the bivariate Skellam distribution, denoted by BSkellam(λ 0 , λ 1 , λ 2 ), if and only if

ε 1 d = U 1 -U 0 , ε 2 d = U 2 -U 0 ,
where U 0 , U 1 and U 2 are three independent random variables such that

U i ∼ Poisson(λ i ) for any i ∈ {0, 1, 2}. By convention, we set U 0 = 0 if λ 0 = 0.
Remark 10 When λ 0 = 0, the distribution of (ε 1 , ε 2 ) is that of a couple of independent Poisson random variables. Now, suppose that λ 0 > 0. Then, the probability mass function of (ε 1 , ε 2 ) is given by

P(ε 1 = x 1 , ε 2 = x 2 ) = e -(λ1+λ2+λ0) λ x1 1 λ x2 2 ∞ i=max(0,-x1,-x2) (λ 0 λ 1 λ 2 ) i (x 1 + i)!(x 2 + i)!i! , (x 1 , x 2 ) ∈ Z 2 .
the mean of (ε 1 , ε 2 ) is (λ 1 -λ 0 , λ 2 -λ 0 ) and the covariance matrix equals

Σ =       λ 1 + λ 0 λ 0 λ 0 λ 2 + λ 0      
.

A detailed study of the bivariate Skellam distribution is postponed to the annex.

Remark 11 Under the above parametric assumption, it follows from Corollary 1 that the B-SINAR(1) process of second kind defined by (4) has a unique stationary solution, and it is also ergodic.

Definition 6 (SBP-SINAR(1) process) A B-SINAR(1) process of second kind with the above parametric assumptions will be denoted by SBP-SINAR(1) (for the first-order Skellam Bivariate Parametric SINAR process) in the following.

Next, two estimation methods are proposed to estimate the parameters of the SBP-SINAR(1) process, i.e., θ = (α 1 , α 2 , λ 0 , λ 1 , λ 2 ).

The method of moments

Based on ( 6), ( 8), and (9) (or (10)), moment estimation is straightforward. Let

X j = 1 n n i=1 X i,j , S 2 j = 1 n -1 n i=1 (X i,j -X j ) 2 , |X| j = 1 n n i=1 |X i,j |, and C 1,2 = 1 n -1 n i=1 (X i,1 -X 1 )(X i,2 -X 2 )
for any j ∈ {1, 2}, and denote the first order sample autocorrelation function by ρXj (1). Then, one obtains

αj = ρXj (1) + 1 2 , ( 13
) λj = λ0 + 2X j (1 -αj ) (14) and λ0 = C 1,2 (1 -ρX1 (1)ρ X2 (1)) ( 15
)
for any j ∈ {1, 2}.

Remark 12 Based on ( 7), one may also introduce other natural estimators for λ j , i.e., λj = (1 -αj )

X j + αj (2S 2 j -|X| j ) (16)
for j ∈ {1, 2}. At first glance, λj seems slightly more complex, and we study the differences between the two estimators by means of a small simulation study in Section 5.

The conditional maximum likelihood estimator

Here, we consider the Conditional Maximum Likelihood (CML) estimator of θ = (α 1 , α 2 , λ 0 , λ 1 , λ 2 ). It is assumed that θ belongs to a compact parameter space Θ. Let x = (x 0 , x 1 , • • • , x n ) be a bivariate sample from SBP-SINAR(1) process with x i = (x i,1 , x i,2 ) for any i ∈ {0, 1, . . . , n}. Thus, the CML estimator is

θCML = argmax θ∈Θ L(x, θ), with L(x, θ) = n t=1 P (X t = x t | X t-1 = x t-1 ) .
Note that the conditional distribution X t | X t-1 in the above equation is equivalent to Γ (•, •), the one-step transition probability of the Markov chain (X t ) t∈Z defined by (4). Hence,

Γ ((a, b), (c, d)) = P (X t = (a, b) | X t-1 = (c, d)) = P (F 1 • c + ε t,1 = a, F 2 • d + ε t,2 = b) .
Let T (j) x be a Binomial(2|x|, α j ) random variable for all x ∈ Z * and any j ∈ {1, 2}. Then, using Remark 8, we distinguish four cases:

-For c = d = 0 holds Γ ((a, b), (0, 0)) = P (ε t,1 = a, ε t,2 = b) ; -For c = 0 and d = 0 holds Γ ((a, b), (c, 0)) = a+|c| k=a-|c| P T (1) c = sign(c)(c + a -k) P (ε t,1 = k, ε t,2 = b) ; -For c = 0 and d = 0 holds Γ ((a, b), (0, d)) = b+|d| k=b-|d| P T (2) d = sign(d)(d + b -k) P (ε t,1 = a, ε t,2 = k) ; -For c = 0 and d = 0 holds Γ ((a, b), (c, d)) = a+|c| k1=a-|c| b+|d| k2=b-|d| P (ε t,1 = k 1 , ε t,2 = k 2 ) × P T (1) c = sign(c)(c + a -k 1 ) P T (2) d = sign(d)(d + b -k 2 ) .
Note that the joint probability mass function of the two innovation process is given in Remark 10.

Remark 13 Indeed, CML estimation is more complicated than the method of moments, but still tractable. Numerical maximization is straightforward with standard statistical packages, but may be computationally costly for long sequences of observations. Remark 14 Since (X t ) t∈Z is a stationary Markov chain, one can obtain asymptotic normality of the parameter estimators from Theorem 2.2 in [START_REF] Billingsley | Statistical Inference for Markov Processes[END_REF] under standard technical assumptions. For more details we refer to [START_REF] Franke | Multivariate first-order integer-valued autoregressions[END_REF].

Simulation study

In this part we shortly investigate some properties of the parameter estimation techniques of the SBP-SINAR(1) process, as presented in Section 3.2. More precisely, we consider a SBP-SINAR(1) model with parameters (λ 0 , λ 1 , λ 2 ) = (7, 3, 8) and (α 1 , α 2 ) = (0.3, 0.8). Figure 1 shows random samples with 200 observations of the two components X 1 and X 2 , the dashed lines correspond to the expected values of the process. The two components are correlated, in this example the estimated correlation coefficient equals 0.206 (p = 0.00332).

The main focus of the Monte Carlo study lies on comparing the properties of the methods of moments and the the maximum likelihood estimator. Therefore, we generated 200 sequences of observations of length n = 50 (respectively 100, 200), using the parameters above. Subsequently, the parameters were estimated by the method of moments and by maximization of the conditional likelihood. With respect to the methods of moments estimator, it should be noted that two methods exist for estimating the parameters λ 1 and λ 2 (see Equations ( 14) and ( 16)). Figure 2 shows Box plots of the estimated parameters. The maximum likelihood estimator is colored dark gray, the methods of moments estimators from Equation ( 14) have white background. Box plots with light gray background result from the estimator of Equation ( 16), and the horizontal dashed lines correspond to the true parameter values. Several aspects are apparent: First, all estimators seem to be consistent and unbiased, only α i , i ∈ {1, 2} might be subject to a slight underestimation for short sequences. Second, the maximum likelihood estimator has the lowest variability. Third, comparing the methods of moments estimators, it seems advisable to prefer Equation ( 16) due to the lower variability of the estimates. Finally, it may be noted that normality of the estimators cannot be rejected for the majority of samples at 5%-level (26 of 36, Shapiro-Wilk).

Analysis of the annual Swedish population rates and the Harvest index data

In the following, we fit different models to Swedish population rates and Harvest index data. The two top panels of Figure 3 shows the time series of the annual Swedish population rates (per thousand population) and the Swedish Harvest index, denoted P t and h t , respectively, for the period 1750 -1849, as reported by [START_REF] Thomas | Social and Economic aspects of Swedish population mouvements[END_REF]. Although such population variation should be real-valued, the recorded rates are rounded values on the one thousand scale. The rates vary from -27 to 16, with a sample mean of 6.69 and sample variance of 34.6. Originally, the Swedish grain harvest was measured on a half-point scale with a total crop failure scored as zero and a superabundant crop scored as nine. In order to avoid observations such as 3.5 we refine the scale of this series by multiplying all the observations by 2. Thus, the Harvest index takes values between 0 and 18, with a sample mean of 10.5 and sample variance of 27.8. The sample ACF and PACF of (P t ) t∈Z and (h t ) t∈Z are given in the middle panels of Figure 3.

Remark 15 Since the annual Swedish rates series has negative observations, the bivariate time series X t = (P t , h t ) τ cannot be fitted by any existing integervalued model based on the classical or generalized thinning operator.

Remark 16 For both series, sample ACF and PACF have values significantly different to zero only at lag 1, if at all. Moreover, the sample correlation equals only 0.0589 (p = 0.56), but the cross-correlation of P t and h t-1 has the value 0.400 (p = 4.11e-05, see also the lower panel of Figure 3). Therefore, one may expect the B-SINAR(1) model to be more adequate for this data series.

Fit of a B-SINAR(1) model

To analyze the present bivariate time series, we consider the following B-SINAR(1) of first kind model

      P t h t       =       F 1,1 • F 1,2 • F 2,1 • F 2,2 •             P t-1 h t-1       +       ε 1 t ε 2 t       . Let θ = (γ 1,1 , γ 1,2 , γ 2,1 , γ 2,2 , µ 1 , µ 2
) τ be the model parameters, where γ i,j and µ j represent the respective means of F i,j and ε j t for any (i, j) ∈ {1, 2} 2 and t ∈ Z. We estimated the parameters using (11). This resulted in γ1,1 = 0.442 Hence, the one-step ahead least squares predictions can be calculated by Pt+1 = 0.442P t + 0.414h t -0.606 and ĥt+1 = -0.119P t + 0.398h t + 6.97.

Remark 18 In general, Pt+1 and ĥt+1 are real-valued, a mapping into the discrete support of the series is obtained by rounding to the nearest integer.

Fit of a SBP-SINAR(1) model

In order to fit the bivariate time series of the annual Swedish population rates and the Swedish Harvest index, we consider the following SBP-SINAR(1) process

      P t h t       =       F 1 • 0 0 F 2 •             P t-1 h t-1       +       ε 1 t ε 2 t      
, with F j as defined in ( 12), and ε t = ε 1 t , ε 1 t follows the bivariate Skellam distribution (see Definition 5), for any j ∈ {1, 2} and t ∈ Z. Let θ = (α 1 , α 2 , λ 0 , λ 1 , λ 2 ) τ be the model parameters, where α 1 and α 2 ) are the parameters of the distribution F 1 (respectively F 2 ), and λ 0 , λ 1 and λ 2 are the parameters of the bivariate Skellam distribution.

We estimated the model parameters by maximization of the likelihood due to the superior performance of the estimator in the previous section. The resulting parameter estimates are α1 = 0.790 .

As before, standard errors, reported on brackets below each estimate, result from a parametric bootstrap with 200 repetitions. Thus, the one-step ahead least squares predictions lead to Pt+1 = (2α 1 -1)P t + λ1 -λ0 = 0.581P t + 2.73, and ht+1 = (2α 2 -1)h t + λ2 -λ0 = 0.195P t + 8.39, Figure 5 shows the observed values of the series of Swedish population rates and Harvest index and the corresponding one-step ahead forecasts resulting from the estimated SBP-SINAR( 1 Summarizing, both proposed models have their small advantages and disadvantages: On the one hand, the B-SINAR(1) of first kind model obtains a better in-sample fit in terms of MAE. On the other hand, the SBP-SINAR(1) constitutes the more parsimonious model, however, the difference in MAE is rather small. This result can be explained by the low correlation but significant cross-correlation of the series at lag 1. A relevant advantage of the SBP-SINAR(1) model may be expected for data with significantly correlation at lag 0, in particular in a simulation context.

Proofs

Proof of Theorem 1

Let • be the Euclidean norm on R 2 , and denote the corresponding operator norm associated to a matrix of size 2 × 2 in the same way. For all x ∈ E = Z 2 , we define the empirical measure µ n by

µ n (•) = 1 n n i=1 π i (x, •),
where π i denotes the i-step transition probability of the Markov chain (X t ) t∈Z defined by (2) for any i ∈ {1, 2, . . . , n}. The idea is to prove that this sequence of probability distribution (µ n ) n∈Z has a sub-sequence converging to some probability distribution µ. By construction, this limit will automatically be a stationary probability distributions of the B-SINAR(1) process. Let V : R 2 → R + be the Lyapunov function defined by V (x) = x 2 . Since ρ(M ) < 1, we have that M k < 1 if k is large enough. For such a k, using (Kachour and Truquet 2011, Lemma 3), we obtain lim sup

x →+∞ E (V (X k ) | X 0 = x) V (x) < 1,
and the Lyapunov criterion is satisfied. Therefore, the sequence (µ n ) n∈N is tight and any cluster point µ is an invariant distribution with µ(V ) < ∞ (see (Duflo 1997, Proposition 2.1.6)). It follows that µ is a measure π-invariant.

From assumption 1, one can deduce that B-SINAR( 1) is an irreducible and aperiodic Markov chain on Z 2 . Therefore, the Markov chain (X t ) t∈Z , as defined in (2), is positive recurrent on E and the invariant measure µ is unique. The existence and uniqueness of a stationary process X = (X 1 , X 2 ) satisfying (2) follows and µ (V ) = E( X 2 ) < ∞. Finally, the ergodicity of the process X easily follows since (X t ) t∈Z is positive recurrent on E.

Proof of Theorem 2

The asymptotic properties of the CLS estimators for the parameters can be established by using the results in [START_REF] Klimko | On conditional least squares estimation for stochastic process[END_REF]. In the following, we present a direct proof inspired by the proof of Theorem 3.2 in [START_REF] Zhang | Inference for INAR(p) processes with signed generalized power series thinning operator[END_REF]. For any (i, j) ∈ {1, 2} 2 , set

M (i,j) n = - 1 2 ∂S n (θ) ∂γ i,j and M (j) n = - 1 2 ∂S n (θ) ∂µ j .
Observe that

√ n θCLS 1 -θ * 1 = 1 n Q -1 1 √ n            M (1,1) n M (1,2) n M (1) n            , (17)
and

√ n θCLS 2 -θ * 2 = 1 n Q -1 1 √ n            M (2,1) n M (2,2) n M (2) n            . ( 18
) Let F n = σ(X 0 , X 1 , • • • , X n ) and R 1 t = X t,1 -γ 1,1 X t-1,1 -γ 1,2 X t-1,2 -µ 1 . Then we have E M (1,1) n | F n-1 = E M (1,1) n-1 + R 1 n X n-1,1 | F n-1 = M (1,1) n-1 + E R 1 n X n-1,1 | F n-1 = M (1,1) n-1 .

Hence, {M

(1,1) n , F n , n ≥ 1} is a martingale. Under the moments assumptions in Theorem 2, we obtain

1 n n t=1 R 1 t X t-1,1 2 a.s. -→ E R 1 t 2 X 2 t-1,1 = Σ 1,1
from the ergodic theorem. Moreover, using (Hall and Heyde 1980, Corollary 3.2), we get 1

√ n M (1,1) n L -→ N (0, Σ 1,1 ) .
Similarly, one can prove that 1

√ n M (1,2) n L -→ N (0, Σ 2,2 ) and 1 √ n M (1) n L -→ N (0, Σ 3,3 ) , Therefore, if follows that 1 √ n c τ            M (1,1) n M (1,2) n M (1) n            L -→ N 0, E (c 1 X t-1,1 + c 2 X t-1,2 + c 3 ) 2 R 1 t 2 for any c = (c 1 , c 2 , c 3 ) τ ∈ (R * ) 3 . From the Cramer-Wold device, one can deduce that 1 √ n            M (1,1) n M (1,2) n M (1) n            L -→ N (0, Σ) . ( 19
)
Thus, in the same way, we find

1 √ n            M (2,1) n M (2,2) n M (2) n            L -→ N (0, Γ ) . ( 20
)
Combining ( 17) and ( 19) on the one hand, and ( 18) and ( 20) on the other hand, we conclude the proof of Theorem 2.

Annex: some complements on the BSkellam(λ 0 , λ 1 , λ 2 ) distribution

Basic properties

First of all, let us recall the definition of the Skellam distribution.

Definition 7 (Skellam(λ 1 , λ 2 ) distribution) Let λ 1 > 0 and λ 2 > 0. We say that the random variable X has the Skellam distribution, denoted by Skellam(λ 1 , λ 2 ), if and only if

X d = U 1 -U 2 ,
where U 1 and U 2 are two independent random variables such that U i ∼ Poisson(λ i ) for any i ∈ {1, 2}.

The probability mass function of X is given by

P(X = x) = e -(λ1+λ2) λ x 1 ∞ i=max(0,-x) (λ 1 λ 2 ) i (x + i)!i! , x ∈ Z.
Further details can be found, e.g., in [START_REF] Skellam | The frequency distribution of the difference between two poisson variates belonging to different populations[END_REF], [START_REF] Karlis | Bayesian modeling of football outcomes: Using the Skellam distribution for the goal difference[END_REF] and [START_REF] Al-Zaid | On the poisson difference distribution inference and applications[END_REF].

Remark 19 Note that X 1 ∼ Skellam(λ 1 , λ 0 ), X 2 ∼ Skellam(λ 2 , λ 0 ) and X 1 - X 2 ∼ Skellam(λ 1 , λ 2 ) if (X 1 , X 2 ) ∼ BSkellam(λ 0 , λ 1 , λ 2 )
. Moreover, X 1 and X 2 are independent if and only if λ 0 = 0.

Lemma 1 below investigates some basic properties of the BSkellam(λ 0 , λ 1 , λ 2 ) distribution.

Lemma 1 Let λ 0 ≥ 0, λ 1 > 0, λ 2 > 0 and (X 1 , X 2 ) ∼ BSkellam(λ 0 , λ 1 , λ 2 ).

Then the characteristic function of (X 1 , X 2 ) is φ (X1,X2) (t 1 , t 2 ) = e λ1(e it 1 -1) e λ2(e it 2 -1) e λ0(e -it 1 -it 2 -1) , (t 1 , t 2 ) ∈ Z 2 .

Remark 20 Using the characteristic function of (X 1 , X 2 ) ∼ BSkellam(λ 0 , λ 1 , λ 2 ), the probability mass function of (X 1 , X 2 ) can be expressed by

P (X 1 = x 1 , X 2 = x 2 ) = 1 (2π) 2 π -π π -π φ (X1,X2) (t 1 , t 2 )e -ix1t1 e -ix2t2 dt 1 dt 2 , (x 1 , x 2 ) ∈ Z 2 .
In particular, this yields the following equality

π -π π -π
e λ1(e it 1 -1) e λ2(e it 2 -1) e λ0(e -it 1 -it 2 -1) e -ixt1 e -iyt2 dt

1 dt 2 = (2π) 2 e -(λ1+λ2+λ0) λ x 1 λ y 2 ∞ i=max(0,-x,-y) (λ 0 λ 1 λ 2 ) i (x + i)!(y + i)!i! , (x, y) ∈ Z 2 .
Lemma 2 below presents a Gaussian approximation result on the BSkellam(λ 0 , λ 1 , λ 2 ) distribution for large λ 1 + λ 0 and λ 2 + λ 0 .

Lemma 2 For λ 1 + λ 0 and λ 2 + λ 0 large enough, we have the approximation

BSkellam(λ 0 , λ 1 , λ 2 ) ≈ N 2             λ 1 -λ 0 λ 2 -λ 0       ,       λ 1 + λ 0 λ 0 λ 0 λ 2 + λ 0            
.

Estimation

In this section, we carry out the estimation of the unknown parameters via the method of moments and the method of the maximum likelihood. Let λ 0 ≥ 0, λ 1 > 0 and λ 2 > 0 be unknown parameters and (X 1,1 , X 2,1 ), . . . , (X 1,n , X 2,n ) be n i.i.d. bivariate random variables with the common distribution BSkellam(λ 0 , λ 1 , λ 2 ).

The method of moments.

Set, for any j ∈ {1, 2},

X j = 1 n n i=1 X j,i , S 2 j = 1 n -1 n i=1 (X j,i -X j ) 2 and C 1,2 = 1 n -1 n i=1 (X 1,i -X 1 )(X 2,i -X 2 ).
Due to Lemma 1, the method of moment gives

-the estimator λj = 1 2 (X j + S 2 j ).
for λ j with j ∈ {1, 2} the estimator λ0 = C 1,2 for λ 0 .

The maximum likelihood method.

For any z = (x, y) ∈ Z 2 and λ = (λ 0 , λ 1 , λ 2 ) ∈ [0, ∞[ 3 , set G(λ, z) = λ x 1 λ y 2 ∞ i=0 (λ 0 λ 1 λ 2 ) i (x + i)!(y + i)!i! .
Lemma 3 For any (x, y) ∈ Z 2 , we have

∂G(λ, z) ∂λ 1 = x λ 1 G(λ, z) + λ 0 λ 1 G(λ, (x + 1, y + 1)), ∂G(λ, z) ∂λ 2 = y λ 2 G(λ, z) + λ 0 λ 2 G(λ, (x + 1, y + 1))
and ∂G(λ, z) ∂λ 0 = G(λ, (x + 1, y + 1)).

For any i ∈ {1, . . . , n}, set

z i = (x 1,i , x 2,i ) ∈ Z 2 , z = (z 1 , . . . , z n ) and λ = (λ 0 , λ 1 , λ 2 ) ∈ [0, ∞[ 3 . The likelihood function is L(λ, z) = n i=1 P(X 1 = x 1,i , X 2 = x 2,i ) = e -n(λ0+λ1+λ2) n i=1 G(λ, z i ).
Therefore the j-th partial derivative of the log-likelihood function is

∂ ln L(λ, z) ∂λ j = -n + n i=1 ∂G(λ,zi) ∂λj G(λ, z i ) .
for any j ∈ {0, 1, 2}. It follows from Lemma 3 that

∂ ln L(λ, z) ∂λ 1 = -n + 1 λ 1 n i=1 x 1,i + λ 0 λ 1 n i=1 G(λ, z i + 1) G(λ, z i ) , ∂ ln L(λ, z) ∂λ 2 = -n + 1 λ 2 n i=1 x 2,i + λ 0 λ 2 n i=1 G(λ, z i + 1) G(λ, z i )
and

∂ ln L(λ, z) ∂λ 0 = -n + n i=1 G(λ, z i + 1) G(λ, z i ) .
Set Z i = (X 1,i , X 2,i ) for any i ∈ {1, . . . , n}. The maximum likelihood estimate λ = ( λ0 , λ1 , λ2 ) satisfies

                       0 = -n + 1 λ1 n i=1 X 1,i + λ0 λ1 n i=1 G( λ, Z i + 1) G( λ, Z i ) , 0 = -n + 1 λ2 n i=1 X 2,i + λ0 λ2 n i=1 G( λ, Z i + 1) G( λ, Z i ) , 0 = -n + n i=1 G( λ, Z i + 1) G( λ, Z i ) . (21)
Therefore

n i=1 G( λ, Z i + 1)/G( λ, Z i ) = n and, a fortiori, λ1 = X 1 + λ0 , λ2 = X 2 + λ0 . (22)
Substituting ( 22) into the last equation in ( 21), one obtains λ0 . Putting λ0 in( 22), we obtain λ1 and λ2 .

Numerical study

In this section we briefly demonstrate some properties of the bivariate Skellam distribution and the parameter estimation procedures. More precisely, we consider two settings (termed Setting 1 and 2 in the following), dealing with a Bivariate Skellam distribution with parameters (λ 0 , λ 1 , λ 2 ) = (4, 2, 3) and (λ 0 , λ 1 , λ 2 ) = (4, 35, 40), respectively. Figure 6 shows random samples with 500 observations each from Setting 1 in the left and Setting 2 in the right panel. The two settings are subject to different correlation, which equals 0.617 (respectively 0.0966), and their estimated values for the samples displayed are 0.628 (p < 2.2e-16) and 0.0975 (p = 0.0292).

In a first Monte Carlo experiment, we investigate the departure from normality in the two settings. For this purpose, we generated 1000 samples for both settings with 1000 observations each and perform the Shapiro-Wilk test for multivariate normality. In Setting 1, having comparably small values of λ i , i = 1, 2, 3, normality is rejected at 5%-level for 87.2% of the samples. However, the corresponding proportion in the second setting reduces to 30.1%, supporting the results of Lemma 2.

In a second Monte Carlo experiment, we investigate the properties of the method of moments and maximum likelihood estimators. To this end, the two parameter sets described above were used for generating 100 sequences of observations of length n = 50, 100, and 200, respectively. Subsequently, the parameters were estimated by the method of moments and by maximization of the likelihood. Figure 7 shows Box plots of the estimated parameters, where the method of moments estimates have white background and the boxes corresponding to the maximum likelihood estimator are coloured gray. The dashed black lines represent the true parameter values. Several aspects are apparent: First, both estimators seem to be consistent. Second, both estimators seem unbiased. Third, the maximum likelihood estimator has a lower variability than the method of moments estimator. Finally, it may be noted that normality of the estimators cannot be rejected for the large majority of samples at 5%-level (30 of 36, Shapiro-Wilk).

Proofs

Proof of Lemma 1. Let (t 1 , t 2 ) ∈ Z 2 . The characteristic function of (X 1 , X 2 ) is then given by

φ (X1,X2) (t 1 , t 2 ) = E(e it1X1+it2X2 ) = E(e it1(U1-U0)+it2(U2-U0) ) = E(e it1U1+it2U2-i(t1+t2)U0 ) = φ U1 (t 1 )φ U2 (t 2 )φ U0 (-t 1 -t 2 )
= e λ1(e it 1 -1) e λ2(e it 2 -1) e λ0(e -it 1 -it 2 -1) .

This ends the proof of Lemma 1.

Proof of Lemma 2.

Set W 1 = X 1 /(λ 1 + λ 0 ) and W 2 = X 2 /(λ 2 + λ 0 ) and let (t 1 , t 2 ) ∈ Z 2 . Using Lemma 1 and e x ∼ 0 1 + x + x 2 /2, we have φ (W1,W2) (t 1 , t 2 ) = φ (X1,X2) (t 1 /(λ 1 + λ 0 ), t 2 /(λ 2 + λ 0 ))
= e λ1(e it 1 /(λ 1 +λ 0 ) -1) e λ2(e it 2 /(λ 2 +λ 0 ) -1) e λ0(e -it 1 /(λ 1 +λ 0 )-it 2 /(λ 2 +λ 0 ) -1)

≈ e λ1(it1/(λ1+λ0)-(1/2)t 2 1 /(λ1+λ0) 2 ) e λ2(it2/(λ2+λ0)-(1/2)t 2 2 /(λ2+λ0) 2 ) × e λ0(-it1/(λ1+λ0)-it2/(λ2+λ0)-(1/2)(t1/(λ1+λ0)+t2/(λ2+λ0)) 2 ) = e it1(λ1-λ0)/(λ1+λ0) e it2(λ2-λ0)/(λ2+λ0) × e -(1/2)(t 2 1 /(λ1+λ0)+t 2 2 /(λ2+λ0)+2t1t2λ0/(λ1+λ0)(λ2+λ0)) .
for any λ 1 + λ 0 and λ 2 + λ 0 large enough. The last term is the characteristic function of a bivariate random variable

Z ∼ N 2             (λ 1 -λ 0 )/(λ 1 + λ 0 ) (λ 2 -λ 0 )/(λ 2 + λ 0 )       ,       1/(λ 1 + λ 0 ) λ 0 /(λ 1 + λ 0 )(λ 2 + λ 0 ) λ 0 /(λ 1 + λ 0 )(λ 2 + λ 0 ) 1/(λ 2 + λ 0 )            
.

Therefore (W 1 , W 2 ) can be approximated by this distribution and, a fortiori,

(X 1 , X 2 ) ≈ N 2             λ 1 -λ 0 λ 2 -λ 0       ,       λ 1 + λ 0 λ 0 λ 0 λ 2 + λ 0             .
Lemma 2 is proved.

Proof of Lemma 3. For any (x, y) ∈ Z 2 , we have

∂G(λ, z) ∂λ 1 = λ x 1 λ y 2 ∞ i=0 (λ 0 λ 1 λ 2 ) i (x + i)!(y + i)!i! = xλ x-1 1 λ y 2 ∞ i=0 (λ 0 λ 1 λ 2 ) i (x + i)!(y + i)!i! + λ x 1 λ y 2 ∞ i=1 i (λ 0 λ 2 ) i λ i-1 1 (x + i)!(y + i)!i! = x λ 1 G(λ, z) + λ x 1 λ y 2 ∞ i=0 (λ 0 λ 2 ) i+1 λ i 1 (x + i + 1)!(y + i + 1)!i! = x λ 1 G(λ, z) + λ 0 λ 1 G(λ, (x + 1, y + 1)).
The other equalities can be proved similarly. Lemma 3 is proved. q q q q q q q q q q q q q q q 0.2 0.3 0.4 0.5 α 1 n = 50 100 200 q q q q q q q q q q q q q q 0.60 0.70 0.80 0.90 α 2 50 100 200 q q q q q 0 5 10 15 λ 0 50 100 200 q q q q q q q q -4 0 2 4 6 8 10 λ 1 n = 50 100 200 M.M. M.M. 2 M.L. q q q q q q q q q q q q 5 10 15 λ 2 50 100 200

Fig. 2 Parameter estimates from 50 sequences of length n = 50, 100, and 200 simulated from a SBB-SINAR(1) model with true parameters (λ 0 , λ 1 , λ 2 ) = (7, 3, 8) and (α 1 , α 2 ) = (0.3, 0.8). Box plots corresponding to the method of moments (M.M.) and maximum likelihood estimator (M.L.) have white and gray background, respectively. The dashed black lines represent the true parameter values. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -5 0 5 -5 0 5

Swedish population rates

x 1

x 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 10 20 30 40 50 20 30 40 50

x 1

x 2 Fig. 6 Two samples with 500 observations each, generated from a bivariate Skellam distribution with true parameters (λ 0 , λ 1 , λ 2 ) = (4, 2, 3) and (λ 0 , λ 1 , λ 2 ) = (4, 35, 40) in the left (respectively right) panel.

q q q q q q q q q q q q 2 3 4 5 6 7 λ 0

Setting 1 n = 50 100 200 q q q q q 0.5 1.5 2.5 3.5 λ 1 50 100 200 q q q q q q q 2 3 4 q q q q q q q q q q 30 35 40

.

  Standard errors, reported on brackets below each estimate, result from a parametric bootstrap with 200 repetitions.Remark 17The spectral radius of the obtained autoregressive matrix, i.e., 1, thus the stationarity condition of Theorem 1 is satisfied.

Figure 4 |

 4 Figure4shows the observed values of the series of Swedish population rates and Harvest index, with the corresponding one-step ahead forecasts based on the B-SINAR(1) model of first kind. The horizontal lines correspond to the sample mean values of both series. The mean absolute error (MAE) for both series equals

|

  Figure5shows the observed values of the series of Swedish population rates and Harvest index and the corresponding one-step ahead forecasts resulting from the estimated SBP-SINAR(1) model. The mean absolute error (MAE) for both series equals MAE(P t ) = 1 99 99

Fig. 1

 1 Figures

Fig. 3 Fig. 4 Fig. 5

 345 Fig. 3 Annual Swedish population rates (left top panel) and annual Swedish Harvest index (right top panel), 1750 -1849, with corresponding ACFs and PACFs (middle panels) and CCF (lower panel)

Fig. 7

 7 Fig.7Parameter estimates from 100 from simulated sequences of length n = 50, 100, and 200. The upper three panels were simulated from a bivariate Skellam with parameters (λ 0 , λ 1 , λ 2 ) = (4, 2, 3), the lower three panels result from the parameters (λ 0 , λ 1 , λ 2 ) = (4, 35, 40) in the left and right panel, respectively. Box plots corresponding to the method of moments (M.M.) and maximum likelihood (M.L.) estimator have white (respectively grey) background. The dashed black lines represent the true parameter values.

140, (3), 667-683.