De cycliciteit van patentaanvragen : alternatieve modellen voor telgegevens : een vergelijking van drie modellen : het Poissonmodel met QML-eigenschap, het negatief binomiale model en het zero-inflated Poissonmodel

Universiteit van Amsterdam Faculteit Economie en Bedrijfskunde

Afstudeerseminar Econometrie & Operationele Research

De cycliciteit van patentaanvragen: alternatieve modellen voor telgegevens Een vergelijking van drie modellen: het Poissonmodel met QML-eigenschap,

het negatief binomiale model en het zero-inflated Poissonmodel

Joosje Goedhart 10738193 30 december 2017

Inhoudelijk begeleider: Dr. J.C.M. van Ophem

Abstract Dit onderzoek vergelijkt drie modellen, het Poissonmodel met QML-eigenschap, het negatief binomiale model en het zero-inflated Poissonmodel, voor onderzoek naar de cy-cliciteit van patentaanvragen en beschrijft een methode om de parameters van deze modellen te schatten met behulp van iteratieve gewogen kleinste kwadaraten. De procycliciteit van patentaanvragen wordt in alle modellen bevestigd. De grote hoeveelheid nullen in de patent-tellingen en het aantal uitschieters met een zeer hoog aantal patenten, maken het lastig een model te vinden dat de data goed beschrijft. Het negatief binomiale model lijkt hiervoor de beste optie, gezien dit model meer rond nul geconcentreerde residuen heeft en gezien de Wald-test die het Poissonmodel verwerpt.

Inhoudsopgave

1 Inleiding 3

2 Theorie en hypotheses 5

2.1 Theorie over de procycliciteit van innovatie . . . 5

2.2 Modellen voor telgegevens . . . 7

2.2.1 Het Poissonmodel met QML-eigenschap voor telgegevens . . . 8

2.2.2 Het negatief binomiale model voor telgegevens . . . 8

2.2.3 Het Zero-inflated Poisson model voor telgegevens . . . 9

3 Empirische methodologie 11 3.1 Iteratieve gewogen kleinste kwadraten . . . 11

3.1.1 Implementatie in R . . . 13

3.1.2 Poissonmodel met QML-eigenschap . . . 13

3.1.3 Negatief binomiaal model . . . 14

3.1.4 Implementatie zero-inflated poisson model (ZIP) . . . 14

3.2 Toetsing van de drie modellen . . . 14

3.2.1 Poissonmodel: Pearson Grootheid . . . 15

3.2.2 Negatief binomiaal model: Wald-test . . . 15

3.2.3 Zero-inflated Poissonomodel: Score-test . . . 16

3.3 Data en variabelen . . . 17

3.3.1 Dataselectie . . . 17

3.3.2 Beschrijvende Statistiek . . . 17

4 Empirische resultaten 19 4.1 Geschatte regressieco¨effici¨enten . . . 19

4.2 Toetsing van de drie modellen . . . 20

5 Conclusie 26

6 Discussie 28

7 Bibliografie 29

1

Inleiding

Schumpeters (1939) voorspelling dat bedrijven investeringen in innovatie concentreren in re-cessieve periodes, waarin de vraag naar hun producten relatief laag is, is nog altijd een on-derwerp van theoretisch en empirisch debat. Recent empirisch onderzoek beargumenteert dat investeringen in Research and Development (R&D) en innovatie juist procyclisch in plaats van anticyclisch zijn, wat inhoudt dat bedrijven investeringen in innovatie strategisch laten samenvallen met periodes van economische groei (Griliches, 1990; Barlevy, 2007; Fabrizio Tsolmon, 2014). Deze discussie is relevant omdat verondersteld wordt dat procyclische uit-gaven aan innovatie hogere kosten tot gevolg hebben in periodes van economische neergang (Barlevy, 2007). Het procyclische patroon van innovatie wordt toegeschreven aan de financi¨ele beperkingen die bedrijven ondervinden in recessieve periodes (Aghion et al., 2012). De hogere sociale zekerheidskosten zouden voor beleidsmakers een reden kunnen zijn om, in recessies, deze financi¨ele beperkingen ten gevolge van overheidsbeleid te verminderen.

In hun onderzoek naar de procycliciteit van innovatie onderzoeken Fabrizio en Tsolmon (2014) de cycliciteit van R&D en van innovatie, gemeten door patentaanvragen. De resultaten van hun onderzoek, waarin aantallen patentaanvragen worden geanalyseerd, suggereren dat zowel R&D als patentaanvragen procyclisch zijn. De verklaring hiervoor zou zijn dat de imple-mentatie van een innovatie imitatie van concurrenten faciliteert, wat de opbrengsten van de originele invoerder reduceert (Fabrizio Tsolmon, 2014). Bedrijven zouden daarom de imple-mentatie of vercommercialisering van innovatie strategisch laten samenvallen met periodes van economische vooruitgang om hun kortetermijnwinst te maximaliseren. Dit procyclische pa-troon van patentaanvragen wordt ook gevonden door Griliches (1990). Daarnaast suggereert zijn onderzoek dat patentaanvragen relatief vroeg in de levensloop van een onderzoeksproject worden gedaan.

ge-bruikten Fabrizio en Tsolmon (2014) een standaard Poissonregressiemodel met quasi- maxi-mumlikelihoodeigenschap (QML) om de relatie tussen het aantal patentaanvragen en veran-deringen in de vraag naar producten te schatten. Dit model veronderstelt dat de variantie van de tellingen evenredig is met het niveau. Er zijn mogelijk twee problemen met dit mo-del: (a) de variantiefunctie is verkeerd gespecificeerd, (b) veel bedrijven patenteren niet elk jaar of zelfs nooit waardoor de patentaanvragen meer nullen kunnen bevatten dan volgens het Poissonmodel met QML-eigenschap. In dit onderzoek worden, naast het Poissonmodel met QML-eigenschap van Fabrizio en Tsolmon (2014), twee modellen vergeleken: (1) een re-gressiemodel met de negatief binomiale verdeling en (2) een zogenaamd zero-inflated Poisson regressie model, waarbij wordt onderzocht welk van deze drie modellen de data het best be-schrijft.

De rest van dit artikel is als volgt ingedeeld. Sectie 2 beschrijft bestaande theorie¨en be-treffende de procycliciteit van investeringen in innovatie, specificeert een aantal alternatieve modellen voor telgegevens en geeft de hypotheses die worden getoetst. In sectie 3 worden de methodologie en de gebruikte dataset beschreven. Paragraaf 4 presenteert en bediscussieert de empirische resultaten. Paragraaf 5 en 6 sluiten af met de conclusies van dit onderzoek, een discussie en suggesties voor verder onderzoek.

2

Theorie en hypotheses

2.1 Theorie over de procycliciteit van innovatie

Barlevy (2007) geeft als theoretische verklaring voor de procycliciteit van R&D-investeringen en innovatie het verlangen van bedrijven om winst te maximaliseren door innovatieve activi-teiten samen te laten vallen met periodes waarin de vraag naar hun producten hoog is. Hij stelt dat de timing van een bedrijf om te investeren in innovatie sterk afhankelijk is van de kortetermijnwinst die de investering tot gevolg heeft, gezien het feit dat anderen, buiten die-genen die de innovatie introduceren, ook profijt van de innovatie kunnen hebben. Dit weegt volgens Barlevy (2007) zwaarder dan de eerder heersende opvatting dat de opportuniteitskos-ten van productiviteitsgroei lager zijn in recessies, wat een aansporing zou zijn om innovatieve activiteiten te ondernemen in recessieve periodes, zoals Schumpeter (1939) veronderstelt.

Het generen van een nieuw product brengt veel verschillende activiteiten met zich mee, waar-onder het investeren in en uitvoeren van R&D en de vercommercialisering van het product, waarmee de implementatie van het nieuwe idee of product wordt bedoeld . Barlevy (2007) neemt aan dat R&D-investeringen en de vercommercialisering van een nieuw product gelijk-tijdig plaatsvinden, omdat bedrijven ongeduldig zouden zijn en niet wachten met het imple-menteren van nieuwe idee¨en. Fabrizio en Tsolmon (2014) wijken af van deze aanname in hun onderzoek naar de procycliciteit van innovatie. Net als Fran¸cois en Lloyd-Ellis (2003) nemen zij aan dat beslissingen betreffende de timing van R&D en innovatie apart van el-kaar worden genomen, waardoor bedrijven in R&D kunnen investeren om innovatieve idee¨en te genereren en de implementatie van deze innovatieve idee¨en vervolgens strategisch kunnen uitstellen. Gezien het ontbreken van data met betrekking tot exacte productintroducties ge-bruiken Fabrizio en Tsolmon (2014) de timing van patentaanvragen van bedrijven als maat voor de vercommercialisering van innovatie. Fabrizio en Tsolmon (2014) stellen dat patent-aanvragen de verdere ontwikkeling en vercommercialisering van een innovatie representeren en de innovatie dichterbij implementatie brengen. Het patentsysteem van de Verenigde Staten

garandeert patentbescherming voor een bepaalde hoeveelheid tijd die start op het moment van de patentaanvraag. Dit geeft bedrijven een drijfveer om de patentaanvraag zo lang mo-gelijk uit te stellen. Deze veronderstelling suggereert dat er een bepaalde vertraging optreedt tussen het investeren in R&D en het aanvragen van een patent. Hall, Griliches en Hausman (1986) onderzoeken de relatie tussen R&D en patenten en trachten de vertragingsstructuur van deze relatie te specificeren. Zij constateren dat er een sterke, bijna gelijktijdige relatie tussen R&D-uitgaven en patentaanvragen bestaat die een aantal kleine vertragingseffecten bevat. Deze vertragingseffecten zijn echter lastig te specificeren en kunnen niet goed worden geschat op basis van de gebruikte data, waardoor zij er niet in slagen deze vertragingsstruc-tuur verder op te helderen. Cohen, Nelson en Walsh (2000) stellen dat de effectiviteit van patenten in het toe-eigenen van de winst van innovatie enorm verschilt tussen industrie¨en. Dit zou kunnen betekenen dat het vertragingseffect tussen het vervullen van een patentaanvraag en investeren in R&D ook sterk verschilt tussen industrie¨en.

Fabrizio en Tsolmon (2014) stellen, net als Barlevy (2007), dat het introduceren van een inno-vatie op de markt imitatie van concurrenten faciliteert, aangezien de innoinno-vatie op dat moment publiekelijk bekend wordt. Om deze reden voorspellen zij dat bedrijven de implementatie van innovatie, gemeten door patentaanvragen, zullen verplaatsen naar periodes van hoge vraag om de omzet van de innovatie te maximaliseren voordat concurrenten hun monopoliepositie wegconcurreren. Daarnaast onderzoeken zij of bedrijven die meer kans hebben op imitatie door concurrenten sterker geprikkeld worden om innovaties te laten samenvallen met periodes van een hoge vraag. Zij concluderen dat het procyclische patroon van innovatie meer uitge-sproken is wanneer de intellectuele eigendomsbescherming van gepatenteerde innovaties, van de industrie waarin het bedrijf opereert, zwakker is.

Fabrizio en Tsolmon (2014) doen hiernaast onderzoek naar de gevoeligheid van de hoeveelheid patentaanvragen voor verandering in de vraag in industrie¨en waarin de verouderingsgraad van producten hoog is. Zij voorspellen dat, in industrie¨en met een hoge mate van

productver-oudering, bedrijven minder snel geneigd zijn de implementatie van innovatie uit te stellen naar periodes van hoge vraag. Deze hypothese, dat patentaanvragen minder sterk procyclisch zijn in industrie¨en met een hoge verouderingsgraad, wordt in het onderzoek van Fabrizio en Tsolmon (2014) echter niet bevestigd.

Aghion et al. (2010) stellen dat de aanwezigheid van kredietbeperkingen ervoor zorgen dat investeringen een procyclisch patroon volgen. Een negatieve schok zou ervoor zorgen dat be-drijven minder interne resources hebben en dat hun toegang tot externe financiering wordt beperkt, waardoor de kans op falen wordt vergroot en investeringen voornamelijk op de korte termijn worden gedaan (Aghion et al., 2010). Fabrizio en Tsolmon (2014) onderzoeken of bedrijven in industrie¨en met een grotere afhankelijkheid van externe financiering een meer procyclisch investeringspatroon hebben, omdat deze bedrijven gevoeliger zouden zijn voor fluctuaties in de beschikbaarheid van externe financiering. Deze hypothese wordt in hun on-derzoek echter niet bevestigd vanwege het ontbreken van significante co¨effici¨enten (Fabrizio en Tsolmon, 2014).

De resultaten van Fabrizio en Tsolmon (2014) tonen een zeer hoge correlatie in de tijd tussen (´e´en-jaar) vertraagde R&D-investeringen en aantallen patenten, wat een bevestiging kan zijn van de onderzochte vertraging tussen R&D-investeringen en patenten van Hall et al. (1986).

2.2 Modellen voor telgegevens

Fabrizio en Tsolmon (2014) toetsen de hypothese dat innovatie, gemeten door aantallen pa-tentaanvragen, procyclisch is, wat betekent dat bedrijven hun patentaanvragen strategisch laten samenvallen met periodes van hoge vraag. Daarnaast onderzoeken zij de invloed van de sterkte van de intellectuele eigendomsbescherming, de hoogte van de productverouderings-graad en de afhankelijkheid van externe financiering in een industrie op de sterkte van deze procycliciteit. De responsvariabele, de aantallen patentaanvragen, is een telling. Winkelmann (2008) stelt dat er overdispersie kan optreden in een Poissonmodel voor tellingen, wat

bete-kent dat de variantie groter is dan het gemiddelde, hetgeen in strijd is met de veronderstelde gelijkheid tussen deze twee parameters in het Poissonmodel. Een Poissonmodel met QML-eigenschap kan een oplossing zijn voor een niet correct gespecificeerde verdeling (Winkelman, 2008).

2.2.1 Het Poissonmodel met QML-eigenschap voor telgegevens

Fabrizio en Tsolmon (2014) gebruiken een Poissonmodel met QML-eigenschap om de respons variabele te modelleren. De belangrijkste veronderstelling in dit model is dat de variantie evenredig is met het niveau, waarbij de evenredigsheidsconstante, ook wel de dispersiepara-meter genoemd, geschat wordt door de Pearson statistische grootheid (Winkelmann, 2008). In het eerste onderdeel van dit onderzoek wordt onderzocht of de Poissonverdeling goed bij de patenttellingen aansluit. Gezien de grote hoeveelheid nullen in de patenttellingen is de verwachting dat deze verdeling de data niet goed beschrijft. Het Poissonmodel met QML-eigenschap wordt geschat door onderstaande likelihoodfunctie te maximaliseren.

Y ∼ Poisson(λ) met λ = exp(X0β) L(y|λ) = log e

−λ_λy

y!
= −λ + ylog(λ)

De standaarddeviaties van de geschatte co¨effici¨enten worden achteraf bepaald met pλφ, met φ de dispersieparameter: φ = 1 n − k X i (yi− ˆλi)2 ˆ λi (2.2.1)

2.2.2 Het negatief binomiale model voor telgegevens

Winkelmann (2008) presenteert een negatief binomiaal model als een model voor tellingen met overdispersie. In dit model is de conditionale variantie, λ, altijd groter dan de conditionele verwachting, λ(1 +_αλ). Dit model kan een alternatief zijn voor het Poisson model wanneer er sprake is van over- of onderdispersie, vanwege de grotere flexibiliteit in het modelleren van de

variantiefunctie. In het tweede onderdeel van dit onderzoek wordt onderzocht of het negatief binomiale model een oplsossing biedt voor mogelijke overdispersie. Hiervoor zal het door Winkelmann (2008) beschreven NegBin II model worden gebruikt, waarbij de verwachting van de negatief binomiale verdeling gespecificeerd wordt in termen van de exogene variabelen. Het NegBin II model is als volgt gedefinieerd:

Y ∼ NegBin(α, λ) met λ = exp(X0β) pdf (y|λ, α) = Γ(α + y) Γ(α)Γ(y + 1)  α λ + α α λ λ + α y

L(y|α, λ) = log(Γ(α + y)) − log(Γ(α)) + αlog( α

α + λ) + ylog( λ α + λ)

(2.2.2)

2.2.3 Het Zero-inflated Poisson model voor telgegevens

Wanneer de data een overmaat aan nullen bevat, dat niet verenigbaar is met een standaard model voor telgegevens, zoals een Poisson of negatief binomiaal model, kunnen zero-inflated modellen een oplossing bieden. Het zero-inflated model onderscheidt twee typen nullen: nullen uit een aditioneel geschatte kans en nullen uit een kansverdeling (Winkelmann, 2008), zoals hieronder te zien. P [Y = y] =      ω + (1 − ω)g(0) als y = 0 (1 − ω)g(y) als y > 0

met g(y) een kansdichtheid voor telgegevens

(2.2.3)

In het zero-inflatedmodel worden twee typen nullen verkregen. In het geval van patentaanvra-gen beschrijft het ene type nullen bedrijven die hebben besloten niet te patenteren, ongeacht of er een uitvinding wordt gedaan, bijvoorbeeld wanneer de patenteffectiviteit laag is. Het andere type nullen representeert bedrijven die patenten registreren, maar dit niet noodzakelijk doen in een periode waarin geen uitvindingen worden gedaan. Lambert (1992) benoemt het eerste type als strategisch en de tweede als toevallig.

van een ZIP-verdeling is gelijk aan λ(1 − ω), en de variantie is gelijk aan λ(1 − ω)(1 + λω). Hieruit volgt dat ook dit model gebruikt kan worden wanneer er sprake is van overdispersie. De hypothese, bij het onderzoek naar het ZIP-model, luidt dat dit model beter bij de data aansluit dan het Poissonmodel met QML-eigenschap van Fabrizio en Tsolmon (2014). De log-likelihoodfunctie die gemaximaliseerd dient te worden is als volgt:

L(γ, ω) =X yi=0 logω + (1 − ω)e−λ+X yi>0 h log(1 − w) − λ + ylog(λ)i

met λ = exp(X0β) en ω de kans op y = 0.

(2.2.4)

De volgende sectie geeft een beschrijving van de data en gaat dieper in op de gebruikte methodologie.

3

Empirische methodologie

In dit onderzoek worden drie modellen vergeleken die de procycliciteit van patentaanvragen onderzoeken. Hierbij wordt de relatie zoals beschreven in Fabrizio en Tsolmon (2014) gebruikt:

E[Pkt|Xit, Zkt−1] = exp[β1RDkt−1+ β2Xit+ β3Mkt−1+ β4Xit× Obs

+ β5Xit× P atef f + β6Xit× ExtF in + τt+ µk]

(3.0.1)

Met Pkthet aantal patenten dat bedrijf k in industrie i in jaar t aanvroeg, Xit de natuurlijke

log van de output in industrie i, Zkt−1 een vector die de natuurlijke log van R&D-uitgaven

van bedrijf k in het jaar t − 1 bevat, Rkt−1 en Mkt−1 een groep van ´e´en jaar vertraagde

be-drijfsniveau controlevariabelen, waaronder de hoeveelheid werknemers (ln EMP), een asset-intensiteit controle variabele (ln PPE) en een variabele voor de omzet (ln SALES). Tot slot zijn µk de bedrijfseffecten en τt de jaareffecten. De kruistermen β4Xit× Obs, β5Xit× P atef f

en β6Xit× ExtF in zijn toegevoegd om de invloed van industriecondities (de

productverou-deringsgraad, patenteffectiviteit en de afhankelijkheid van externe financiering) op de gevoe-ligheid van patentaanvragen op een verandering in de vraag te meten. Bij de analyse van model (3.0.1) wordt gecorrigeerd voor verschillen tussen jaren en bedrijven, dat wil zeggen dat er vaste jaar- en bedrijfseffecten worden meegenomen. Dit model wordt op drie manieren geschat, waarbij de responsvariabele wordt gemodelleerd met (1) de Poissonverdeling, (2) de negatief binomiale verdeling en (3) de zero-inflated Poissonverdeling. Onderzoek naar de rela-tie tussen de gestandaardiseerde residuen, de residuen gedeeld door de standaard deviarela-ties, en de gefitte waarden van deze modellen moeten in combinatie met een aantal parametertoetsen van de modellen aangeven welke verdeling de data het best representeert.

3.1 Iteratieve gewogen kleinste kwadraten

Om vergelijking (3.0.1) te schatten met het Poisson-, het negatief binomiale en het ZIP-model dienen de log-likelihoodfuncties van deze ZIP-modellen gemaximaliseerd te worden. Bij

deze schatting dient er gecorrigeerd te worden voor verschillen tussen bedrijven en jaren. De hoeveelheid bedrijven is echter dusdanig groot dat Stata en R geen standaard pakket-ten bieden die een maximum likelihood optimalisatie kunnen uitvoeren voor de beschreven modellen, met uitzondering van het Poissonmodel met QML-eigenschap in Stata. Om deze reden wordt de optimalisatie van de log-likelihoodfuncties gedaan door Newton-Rhapsons methode voor gegeneraliseerde lineaire modellen te herschrijven als een iteratieve gewogen kleinste kwadratenmethode, zoals voor het eerst beschreven in Nelder en Wedderburg (1972). Newton-Rhapson maximaliseert een log-likelihoodfunctie F (β) door een initi¨ele schatting β te verbeteren met behulp van β1 = β0− H−1g met H de Hessiaan en g de gradi¨ent van F (β)

ge¨evalueerd in β0. Dit iteratieschema kan voor gegeneraliseerde lineaire modellen (GLM)

worden omgeschreven naar een iteratief gewogen kleinste kwadratenmethode, ook wel bekend als Fishers scoremethode. (McCulllagh en Nelder, 1989). Voor een GLM met L(µi) de

log-likelihood voor observatie i, µi = h(ηi) de inverse van de linkfunctie die µi aan de lineaire

voorspeller ηi relateert, en met ηi = β0+ β1xi1+ ... + βkxik kunnen de Hessiaan en gradi¨ent

als volgt bepaald worden.

∂L ∂βj =X i ∂L(µi) ∂βj =X i ∂L(µi) ∂µi ∂µi ∂ηi ∂ηi βj =X i ∂L(µi) ∂βj h0(ηi)xij = X i sixij = Xj0s

met si= L0(µi)h0(ηi) en Xj de j-de kolom van de design matrix X. De vector van

eerste orde afgeleiden wordt dan gegeven door X’s, met s de scorevector met elementen si.

Op vergelijkbare wijze kan de Hessiaan verkregen worden als: H = −X0W X met W een diagonaalmatrix met gewicht wi voor observatie i, en

wi = −L00(µi)[h0(ηi)]2− L0(µi)h00(ηi)

(3.1.1)

Met behulp van (3.1.1) kan de Newton-Rhapsonbenadering worden geschreven als β1= β0+

uitdrukking worden afgeleid.

(X0W X)β1= X0W yadjusted (3.1.2)

Uitdrukking (3.1.2) is gelijk aan de gewogen normaalvergelijking van een lineair regressiemo-del. Hieruit volgt dat β1 verkregen kan worden door een regressie van yadjustedop de covariaten

X met gewichten w. Dit proces kan iteratief herhaald worden totdat convergentie optreedt (McCullagh en Nelder, 1989).

3.1.1 Implementatie in R

Voor de drie modellen wordt Fishers scoremethode ge¨ımplementeerd in R. Vergelijking (3.1.2) dient geschat te worden met vaste bedrijf- en jaareffecten. Dit wordt gedaan met de R-functie felm, fixed effects linear models, die gebruikt kan worden om lineaire modellen met meerdere vaste groepeffecten te schatten. Hierbij wordt gebruik gemaakt van de methode van ’alternating projections’ om de groepeffecten uit de normaalvergelijking te vegen voordat de overgebleven co¨effici¨enten worden geschat met OLS. Wanneer aan deze functie gewichten worden toegekend wordt de gewogen kleinste kwadratenmethode uitgevoerd. De volledige R-scripts zijn te vinden in de appendix. Daarnaast is voor elk model een specificatie van de correctheid van de iteratief gewogen kleinste kwadratenmethode gegeven, door deze methode, voor een kleine dataset, te vergelijken met de standaard R of Stata methodes.

3.1.2 Poissonmodel met QML-eigenschap

Voor de implementatie van het Poissonmodel met QML-eigenschap, zoals in Fabrizio en Tsol-mon (2014), in R wordt de log-likelihoodfunctie (2.2.10 gebruikt. De lineaire regressie (3.1.2) is de regressie zoals beschreven in (3.0.1) met vaste bedrijf- en jaareffecten en gewichten w, geschat met de R-functie felm. De lineaire parameters β worden geschat met behulp van Fis-hers scoremethode, zoals beschreven in (3.1.2). De afleiding van de Hessiaan en de grad¨ıent van het Poissonmodel met QML-eigenschap zijn te vinden in Appendix A en de vergelijking

van de iteratief gewogen kleinste kwadratenmethode voor dit model met standaard R en Stata methodes in Appendix E en F.

3.1.3 Negatief binomiaal model

Bij het implementeren van het NegBin II model wordt log-likelihoodfunctie (2.2.2) gebruikt. De lineaire regressiefunctie (3.1.2) is dezelfde als in het Poissonmodel. In het NegBin II-model dienen, naast de lineaire parameters β, ook de α-parameter te worden geschat. Dit wordt gedaan met behulp van een zogenaamd twee-stapsalgoritme waarin de volgende stappen worden herhaald: (1) gegeven een schatting van α, schat de lineaire parameters uit λ = exp(X0β) met 3.1.2., en (2) gegeven de lineaire parameters, schat de parameter α. De afleiding van de Hessiaan en gradi¨ent gespecificeerd zoals in 3.1.1 zijn te vinden in Appendix B. Een vergelijking van de iteratief gewogen kleinste kwadratenmethode met standaard R-methodiek is voor het negatief binomiale model te vinden in Appendix H.

3.1.4 Implementatie zero-inflated poisson model (ZIP)

Het ZIP-model wordt op dezelfde manier als het NegBin II ge¨ımplementeerd. De log-likelihood-functie (2.2.4) wordt gemaximaliseerd door iteratief gewogen kleinste kwadraten toe te passen, waarbij dezelfde lineaire regressiefunctie wordt gebruikt als bij het Poissonmodel. De parame-ters die hier geschat worden zijn de lineaire parameparame-ters uit λ = exp(X0β) en ω, de additionele kans op een telling gelijk aan nul, waarvoor wederom het twee-stapsalgoritme wordt gebruikt. De afleiding van de Hessiaan en de gradi¨ent van de log-likelihoodfunctie van het ZIP-model zijn te vinden in Appendix C. Een vergelijking van de iteratief gewogen kleinste kwadraten-methode met standaard R-methodiek is voor het ZIP-model te vinden in Appendix J.

3.2 Toetsing van de drie modellen

Om veronderstellingen te doen over welk van de drie onderzochte modellen, zoals beschreven in sectie 2.2, de data het best beschrijft, worden de modellen getoetst. Voor het Poissonmodel met QML-eigenschap en het negatief binomiale model wordt getoetst of er sprake is van

overdispersie, dat wil zeggen of de geschatte variantie van het model groter is dan de geschatte verwachting. Voor het zero-inflated Poissonmodel wordt getoetst of de geschatte kans op nullen significant van nul verschilt. In de rest van dit onderdeel wordt besproken welke toetsen worden gebruikt om de parameters van de drie modellen te testen.

3.2.1 Poissonmodel: Pearson Grootheid

Winkelmann (2008) beschrijft de Pearson Grootheid als optie om de goodness-of-fit van het standaard Poissonmodel voor tellingen te meten, waarbij P in (3.2.1) de Pearson grootheid voorstelt. Winkelmann (2008) stelt dat P 6= (n − k) een indicatie is van een misspecificatie van de verwachting van het Poissonmodel. Wanneer dit het geval is, kan er sprake van onder-of overdispersie en kan een Poissonmodel met QML-eigenschap een oplossing bieden, waarbij φ in 2.2.1 gelijk is aan de Pearson grootheid gedeeld door n-k.

P = n X i=1 (yi− ˆλi)2 ˆ λi

, waarbij y het standaard Poissonmodel voor tellingen volgt, ˆ

λ = exp(X0β) en met n het aantal waarnemingenˆ

(3.2.1)

3.2.2 Negatief binomiaal model: Wald-test

Om het Poissonmodel tegen het negatief binomiale model te toetsen met nulhypothese H0 :

θ = 0 versus H1 : θ 6= 0, met θ = _α1 uit het negatief binomiale model in 2.2.2, beschrijft

Winkelmann (2008) een aantal mogelijke toetsen, waaronder de Wald-test. Onder de nulhy-pothese volgt de Wald grootheid, zoals hieronder gespecificeerd, een chi-kwadraatverdeling met 1 vrijheidsgraad (Winkelmann, 2008). Winkelmann merkt op dat de Wald-toets gebruikt kan worden wanneer maximumlikelihood als schattingsmethode wordt gebruikt. Aangezien iteratief gewogen kleinste kwadraten exact dezelfde resultaten oplevert als maximumlikeli-hood, zoals te zien in Appendix H, kan deze toets ook gebruikt worden wanneer de schatting wordt gedaan aan de hand van iteratief gewogen kleinste kwadraten. De Wald grootheid is

als volgt gedefinieerd:

W = (ˆθ − θ)

2

V ar(ˆθ) waarbij onder de nulhypothes geldt W ∼ χ

2_{(1). (3.2.2)}

Omdat in dit onderzoek niet θ, maar α = 1_θ wordt geschat, luidt de nulhypothese hier H0 :

α = 0 versus H1 : α 6= 0. Om deze hypotheses te toetsen, wordt de variantie uit de Wald

grootheid omgeschreven met behulp van een eerste orde Taylorbenadering, zoals hieronder beschreven.

f (x) ≈ f (µ) + (x − µ)f0(µ) ⇒ V ar[f (x)] ≈ V ar(f (µ) + (x − µ)f0(µ)) = V ar(x)[f0(µ)]2 Dus er geldt: V ar(1

θ) ≈ V ar(α) 1

α4, waarmee de Wald-grootheid uit 3.2.2 omgeschreven kan

worden naar: W = (α − α0)

2

V ar(α)_α14

, waarbij W ∼ χ(1) onder de nulhypothese .

(3.2.3)

3.2.3 Zero-inflated Poissonomodel: Score-test

Om te toetsen of de geschatte kans op nullen in het Zero-inflated Poissonmodel significant van nul verschilt, kan de Scoretest gebruikt worden (Mulhally, 1986), waarbij de hypotheses H0 : ω = 0 versus H1 : ω 6= 0 worden getoetst. De score grootheid ξ, zoals hieronder

gespecificeerd, volgt een asymptotische chi-kwadraatverdeling met ´e´en vrijheidsgraad onder de nulhypothese (Mulhally, 1986). De score grootheid wordt berekend met behulp van de scorevector en de informatiematrix uit het gerestricteerde model, met ω = 0.

ξ = s( ˆξ)0Υ( ˆξ)−1s( ˆξ), met s( ˆξ) de scorevector en Υ( ˆξ) de informatiematrix, beide ge¨evalueerd in de ML-schattingen van het gerestricteerde model.

(3.2.4)

In het ZIP-model levert iteratief gewogen kleinste kwadraten dezelfde co¨effici¨entwaardes op als wanneer een ML-schatting wordt gedaan. Hierbij moet worden opgemerkt dat de Hessiaan

van de ZIP-loglikelihoodfunctie niet blokdiagonaal is, waarmee bij iteratief gewogen kleinste kwadraten geen rekening wordt gehouden. In Appendix J is te zien dat de Hessiaan in werkelijkheid groter is dan op basis van iteratief gewogen kwadraten wordt verondersteld. Echter valt het opstellen van de Hessiaan inclusief niet-diagonale termen buiten de scope van dit onderzoek en wordt de Scoretest uitgevoerd met de blokdiagonale Hessiaan uit de iteratief gewogen kleinste kwadratenmethode.

3.3 Data en variabelen 3.3.1 Dataselectie

In dit onderzoek wordt de paneldataset gebruikt zoals gegeven in Fabrizio en Tsolmon (2014). Voor een uitgebreide specificatie van de constructie van deze dataset wordt naar Fabrizio en Tsolmon (2014) verwezen. De gebruikte gegevens beslaan de periode 1975 tot 2002 en zijn afkomstig uit de Verenigde Staten. Belangrijk om op te merken is dat Fabrizio en Tsolmon (2014) in hun analyses geen bedrijven meenemen die geen enkele keer patenteren in alle jaren van de steekproef. In dit onderzoek worden deze waarnemingen wel meegenomen. Bedrijven die over de hele steekproef niet investeren in R&D worden niet meegenomen in onderstaande analyses, omdat het zeer onwaarschijnlijk is dat bedrijven die nooit investeren in R&D pa-tentaanvragen zullen doen. Na het verwijderen van observaties met missende waarden bevat de dataset in totaal 49.1894 bedrijf-jaarobservaties, die 705795 patenten door 4444 bedrijven in 100 verschillende industri¨en representeren.

3.3.2 Beschrijvende Statistiek

De beschrijvende statistiek van de variabelen gebruikt in de analyses wordt in figuur 1 weerge-ven. Gemiddeld vragen bedrijven 14 patenten per jaar aan, terwijl de mediaan van patenten gelijk is aan 1. De patentwaarnemingen bestaan voor 49,19% uit nullen en bevatten een aantal hoge uitschieters zoals te zien in figuur

Statistic Mean St. Dev. Min Median Max NumPats 14.412 76.213 0 1 2,655 ln R&D (lag) 5.195 3.038 0.000 5.606 13.761 ln R&D ( 2lag) 5.007 3.089 0.000 5.437 13.761 ln Sales (lag) 8.775 3.273 0.000 9.051 16.841 ln PPE (lag) 7.380 3.083 0.000 7.420 16.077 ln Emp (lag) 0.995 1.221 0.000 0.445 6.777 ln Output 9.011 1.086 5.790 8.855 12.368 Patenteffectiviteit 0.462 0.145 0.050 0.450 0.821 Verouderingsgraad 2.870 0.345 1.564 2.867 4.000 External Finance 0.644 1.632 −2.934 −0.013 4.762 year 1,975 2,002

Tabel 1: Beschrijvende Statistiek

4

Empirische resultaten

Figuur 4 rapporteert de resultaten van de schattingen van de patentvergelijking (3.0.1) op basis van drie verschillende modellen. In de schatting wordt gecorrigeerd voor verschillen tussen bedrijven en jaren, dat wil zeggen dat vaste bedrijf- en jaareffecten zijn opgenomen in het model. Het Poissonmodel met QML-eigenschap schat de relaties op dezelfde manier als Fabrizio en Tsolmon (2014) en is opgenomen als referentiekader. De schattingen uit het negatief binomiale model (NegBin II) en het zero-inflated Poissonmodel (ZIP) zijn uitgevoerd zoals beschreven in sectie 2 en 3.

4.1 Geschatte regressieco¨effici¨enten

Uit tabel 3 en 4 blijkt dat de geschatte effecten, hier de richting van de co¨effici¨enten, niet verschillen over de modellen. In tabel 3 is te zien dat de co¨effici¨ent die de invloed van de industrieoutput op de hoeveelheid patentaanvragen schat, sterk significant en positief is in de drie modellen. Dit bevestigt dat innovatie, gemeten door patentaanvragen, procyclisch is voor de onderzochte steekproef, wat in overeenstemming is met de waarde zoals geschat in Fabrizio en Tsolmon (2014). Daarnaast is te zien dat de geschatte standaarddeviaties van deze co¨effici¨ent een factor 5 of 10 lager zijn voor de drie modellen in tabel 3.

De hypothese dat patentaanvragen van bedrijven in industri¨en met een relatief hoge paten-tefficitiviteit minder gevoelig zijn voor fluctuatie in de conjunctuurbeweging wordt bevestigd door de negatieve waarde van de significante interactieco¨effici¨ent (ln Output x Patenteffec-tiviteit) in tabel 4. Voor het Poissonmodel met QML-eigenschap en het ZIP-model komt deze schatting sterk overeen met de schatting van Fabrizio en Tsolmon (2014), maar voor het negatief binomiale model is deze co¨effici¨ent in vergelijking met Fabrizio en Tsolmon (2014) lager geschat. Opnieuw zijn de standaard deviaties van deze co¨effici¨ent voor de drie modellen lager dan in Fabrizio en Tsolmon (2014).

De geschatte co¨effici¨ent van ´e´en-jaar vertraagde R&D-investeringen is in zowel tabel 3 als 4 positief en sterk significant, wijzend op een hoge correlatie over de tijd tussen R&D en

patentaanvragen.

In tabel 4 is te zien dat de term die de interactie meet tussen de verouderingsgraad en output (ln Output x Verouderingsgraad) in het Poissonmodel met QML-eigenschap en het NegBin II-model niet significant is, wat ook zo is in het onderzoek van Fabrizio en Tsolmon (2014). De negatieve waarde van deze significant geschatte term in het ZIP-model suggereert kleinere fluctuaties in patentaanvragen bij veranderingen in de industrie-output in industri¨en met een hogere verouderingsgraad.

De positieve waarde van de sterk significante schatting van de interactie term (ln Output x External Finance) bevestigt dat bedrijven in industrie¨en met meer mogelijkheid tot externe financi¨ering gevoeliger zijn voor fluctuaties in de vraag. In het onderzoek van Fabrizio en Tsolmon (2014) wordt deze co¨effici¨ent ook positief geschat, maar zijn de schattingen niet significant.

Concluderend kan over tabel 4 gezegd worden dat deze resultaten empirisch bewijs leveren voor de veronderstelling dat patentaanvragen sterker procyclisch zijn in industri¨en met hoge externe financiering en minder procyclisch in industri¨en met een hoge intellectuele eigendom-bescherming.

4.2 Toetsing van de drie modellen

In tabel 2 zijn de geschatte parameterwaarden van de drie modellen, zoals besproken in sectie 2.2, en de toetsgrootheden, zoals besproken in sectie 3.2, weergeven. In deze tabel is te zien dat de Pearson grootheid van het standaard Poissonmodel voor tellingen gelijk is aan 249425.6, terwijl n-k gelijk is aan 44494. Dit impliceert, zoals beschreven in (3.2.1), dat het standaard Poissonmodel een niet-correct gespecificeerde verwachtingsfunctie heeft en dat er overdispersie optreedt. De dispersieparameter van het Poissonmodel met QML-eigenschap wordt vervolgens geschat op 5.6, wat impliceert dat de variantie van het gebruikte model op ruim vijf keer zo groot als de verwachting wordt geschat. In tabel 2 is te zien dat de α-co¨effici¨ent van het negatief binomiale model geschat wordt op 2.28 en dat de nulhypothese van de Wald-toets verworpen wordt. Dit betekent dat er significant bewijs is dat α ongelijk

is aan nul en dat het negatief binomiale model over het Poissonmodel wordt verkozen. De additieve kans op nullen ω wordt door het zero-inflated Poissonmodel geschat op 0.09, zoals te zien in tabel 2. Uit de waarde van de score-toets blijkt echter dat er geen significant bewijs is dat ω verschilt van nul. Bij het verkrijgen van deze toetsgrootheid is gebruik gemaakt van een Hessiaan die kleiner is dan de werkelijke Hessiaan, waardoor de resultaten in tabel 12 geen werkelijk bewijs leveren voor de veronderstelling dat ω niet significant van nul verschilt.

4.3 Vergelijking van de modellen

De geschatte co¨effici¨entwaardes hebben dezelfde richting over de drie modellen, waardoor dezelfde veronderstelilngen kunnen worden gemaakt op basis van de drie modellen. Deze ver-onderstellingen zijn in overeenstemming met de resultaten van Fabrizio en Tsolmon (2014), behalve dat de hier onderzochte modellen kleinere standaarddeviaties hebben. Dit kan ener-zijds te wijden zijn aan de andere manier van dataopschoning die is toegepast en anderener-zijds aan de twee alternatieve modellen die zijn onderzocht: het negatief binomiale en het ZIP-model. In figuur 7 worden voor elk model de Pearson gestandaardiseerde residuen geplot tegen de gefitte patentaanvragen, waarbij de standaardisatie gegeven wordt in (4.3.1). Voor het Poissonmodel met QML-eigenschap en het zero-inflated Poissonmodel, respectievelijk fi-guur 7(A) en 7(C), is te zien dat de gestandaardiseerde residuen zeer hoge waarden aannemen wanneer de gefitte waarde rond nul ligt. Dit kan impliceren dat deze modellen de nullen in de patenttellingen niet goed representeren. Uit figuur 7(B) blijkt dat de gestandaardiseerde residuen van het negatief binomiale model meer rondom nul gecentreerd zijn, waardoor dit model over het algemeen beter past. Hier is ook te zien dat de residuen voor de uitschieters, d.w.z. de zeer hoge tellingen, ongeveer gelijk zijn aan nul, wat mogelijk impliceert dat de parameterschattingen van het negatief binomiale model in hoge mate worden bepaald door de hoge waarnemingen. Zo is te zien dat dit model gefitte patentaanvragen heeft van rond de 4000, terwijl de maximale hoeveelheid werkelijke patentaanvragen gelijk is aan 2655. Het ZIP-model neemt hetzelfde patroon aan als het negatief binomiale model en toont dezelfde gevoeligheid voor waarnemingen met een hoge waarde. Het patroon van figuur 7(B) en 7(C)

duidt op niet goed gespecificeerde variantiefuncties voor het negatief binomiale en het zero-inflated Poissonmodel. De variantiefunctie van het Poissonmodel met QML-eigenschap lijkt beter gespecificeerd te zijn, in de zin dat de residuen min of meer homogeen zijn. Echter heeft dit model zeer grote residuen in vergelijking met het negatief binomiale model.

Quasipoisson: ei,std=

yi− ˆλi

q ˆ φ ˆλi

met ˆφ de geschatte dispersieparameter, geschat op 5.6.

Negatief Binomiaal: ei,std=

yi− ˆλi q ˆ λi(1 + ˆ λi ˆ α) met ˆα geschat op 2.28.

Zero-inflated Poisson: ei,std=

(1 − ˆω)(yi− ˆλi)

q ˆ

λi(1 − ˆω)(1 + ˆω ˆλi)

met ˆω de geschatte kans op

y=0, geschat op 0.09.

Model Parameter Toetsgrootheid Toetswaarde onder H0

Poisson φ = 5.6ˆ Pearson = 249426 (n-k) = 44494 Negatief Binomiaal α = 2.280324ˆ Wald = 5265 χ2_0.05(1) = 3.841 Zero-inflated Poisson ω = 0.0903ˆ Score = 0.00355 χ2_0.05(1) = 3.841

Tabel 2: : Geschatte parameterwaardes en toetsgrootheden zoals besproken in sectie 3.2.

Afhankelijke variabele: NumPats QuasiPoisson NegBin II ZIP

(1) (2) (3) ln R&D (lag) 0.031∗∗∗ 0.052∗∗∗ 0.052∗∗∗ (0.002) (0.006) (0.002) ln R&D (2lag) 0.020∗∗∗ 0.007 0.007∗∗∗ (0.002) (0.005) (0.002) ln Sales (lag) 0.175∗∗∗ -0.036∗∗∗ -0.036∗∗∗ (0.004) (0.008) (0.005) ln Emp (lag) 0.156∗∗∗ 0.317∗∗∗ 0.317∗∗∗ (0.004) (0.015) (0.004) ln PPE (lag) 0.134∗∗∗ 0.133∗∗∗ 0.133∗∗∗ (0.004) (0.010) (0.004) ln Output 0.238∗∗∗ 0.218∗∗∗ 0.218∗∗∗ (0.002) (0.013) (0.002) Significiantieniveaus: ∗p<0.1;∗∗p<0.05; ∗∗∗p<0.01

Tabel 3: Geschatte co¨efficienten van vergelijking (3.0.1), zonder kruistermen met patenteffec-tiviteit, de verouderingsgraad en externe financiering, waarbij de log-likelihoodfunctie van elk model geoptimaliseerd is met behulp van iteratief gewogen kleinste kwadraten.

Afhankelijke variabele: NumPats QuasiPoisson NegBin II ZIP

(1) (2) (3) ln R&D (lag) 0.029∗∗∗ 0.052∗∗∗ 0.026∗∗∗ (0.005) (0.006) (0.002) ln R&D (2lag) 0.021∗∗∗ 0.006 0.015∗∗∗ (0.004) (0.005) (0.002) ln Sales (lag) 0.185∗∗∗ -0.038∗∗∗ 0.189∗∗∗ (0.011) (0.008) (0.005) ln Emp (lag) 0.196∗∗∗ 0.321∗∗∗ 0.198∗∗∗ (0.009) (0.015) (0.004) ln PPE (lag) 0.088∗∗∗ 0.132∗∗∗ 0.067∗∗∗ (0.011) (0.010) (0.005) ln Output 0.836∗∗∗ 0.487∗∗∗ 0.844∗∗∗ (0.044) (0.102) (0.019) ln Output x Patenteffectiviteit -1.482∗∗∗ -0.512∗∗∗ −1.464∗∗∗ (0.091) (0.145) (0.039) ln Output x Verouderingsgraad -0.019 -0.040 −0.021∗∗∗ (0.018) (0.035) (0.008) ln Output x External Finance 0.094∗∗∗ 0.135∗∗∗ 0.076∗∗∗

(0.010) (0.019) (0.004)

Significiantieniveaus: ∗p<0.1; ∗∗p<0.05; ∗∗∗p<0.01

Tabel 4: Geschatte co¨efficienten van vergelijking (3.0.1), waarbij de log-likelihoodfunctie van elk model geoptimaliseerd is met behulp van iteratief gewogen kleinste kwadraten.

Figuur 4: Gestandadiseerde residuen geplot tegen gefitte waardes uit schattingen door maxi-malisatie van de log-likelihoodfuncties van de volgende modellen: (A) Poissonmodel met QML-eigenschap, (B) Negatief Binomiaal model, (C) Zero-inflated Poissonmodel.

5

Conclusie

Dit onderzoek vergelijkt drie modellen voor onderzoek naar de cycliciteit van patentaanvra-gen en beschrijft een methode om de parameters van deze modellen te schatten met behulp van iteratieve gewogen kleinste kwadraten (McCullagh & Nelder, 1989). Voor het negatief binomiale model en het zero-inflated Poissonmodel is onderzocht of deze modellen de patent-tellingen beter beschrijven dan het Poissonmodel met QML-eigenschap zoals ge¨ımplementeerd in Fabrizio en Tsolmon (2014) en of deze modellen het mogelijke overschot aan nullen in de patentaanvragen goed representeren.

Uit de Pearson grootheid (Winkelmann, 2008) van het geschatte Poissonmodel blijkt dat de verwachtingsfunctie van het standaard Poissonmodel niet goed past. Het resultaat van de Wald-test, die het Poissonmodel tegen het negatief binomiale model toetst, verwerpt de nulhypothese dat de α-parameter van het negatief binomiale model niet significant van nul verschilt, waardoor verondersteld kan worden dat het negatief binomiale model over het stan-daard Poissonmodel wordt geprefereerd. De Score-test in het zero-inflated Poissonmodel toont aan dat de geschatte kans op aditieve nullen niet significant van nul verschilt, waardoor geen bewijs wordt gevonden voor een daadwerkelijk overschot aan nullen in de data.

Uit de resultaten blijkt dat de richting van de geschatte co¨effici¨enten weinig verschilt tus-sen de drie modellen en in overeenstemming is met de richting zoals gevonden in Fabrizio en Tsolmon (2014). De resultaten zijn consistent met het verwachte procyclische patroon van innovatie, gemeten door patentaanvragen. Daarnaast wordt aangetoond dat patentaanvragen sterker procyclisch zijn in industri¨en met lage patenteffectiviteit in vergelijking met industri¨en met hoge patenteffectiviteit. Dit resultaat is consistent met de economische theorie die stelt dat bedrijven innovatie strategisch laten samen vallen met periodes van een hoge vraag, om de korte termijn omzet te maximaliseren (Barlevy, 2007; Fabrizio en Tsolmon, 2014). Wanneer patenten effectiever zijn neigen bedrijven minder naar het matchen van innovatie met periodes

van hoge vraag, vanwege het verminderde risico op immitatie van concurrenten (Fabrizio en Tsolmon, 2014). Tot slot wordt aangetoond dat patentaanvragen sterker procyclisch zijn in industrie¨en die sterker afhankelijk zijn van externe financiering. Dit resultaat is in overeen-stemming met het onderzoek van Aghion et al. (2010) waarin wordt gesteld dat een negatieve shock de kredietbeperkingen in een industrie vergroot en dat bedrijven die sterk afhankelijk zijn van externe financiering vooral kortetermijn innovaties zullen uitvoeren.

Het Poissonmodel met QML-eigenschap geeft over de gehele linie een overmaat aan slecht gefitte waardes. Dit lijkt te impliceren dat de verdelingsfunctie van dit model niet goed bij de data past. Het zero-inflated Poissonmodel geeft minder grote residuen. Echter lijkt het gebruik van dit model overbodig, gezien de geschatte kans op aditieve nullen die niet signi-ficant van nul verschilt. Daarnaast lijkt het onaannemelijk dat de grote waarnemingen een Poissonverdeling zonder overdispersie volgen. De gestandaardiseerde residuen van het negatief binomiale model zijn het meest rond nul gecentreerd, maar de residuen zijn niet homogeen verdeeld.

De resultaten van dit onderzoek leveren empirisch bewijs voor de procycliciteit van inno-vatie, gemeten door patentaanvragen. De combinatie van de grote hoeveelheid nullen in de patenttellingen en een aantal hoge uitschieters maken het lastig een model te vinden dat de data goed beschrijft. Het negatief binomiale model lijkt hiervoor de beste, of minst slechte, optie gezien de relatief meer rond nul geconcentreerde residuen die dit model heeft en ge-zien het resultaat van de Wald-test die het Poissonmodel tegen het negatief binomiale model toetst.

6

Discussie

Dit onderzoek kent twee belangrijke kanttekeningen die hieronder worden besproken. Aller-eerst bevat de data een kleine hoeveelheid waarnemingen met een zeer hoog aantal patentaan-vragen. Deze groep uitschieters draagt deels bij aan de complexiteit van het vinden van een model dat de data goed beschrijft. Vervolgonderzoek zou zich kunnen richten op het vinden van een goed passend model, waarbij deze uitschieters niet worden meegenomen in de analyse. Ten tweede worden er in dit onderzoek en in Fabrizio en Tsolmon (2014) conclusies getrokken over de cycliciteit van innovatie, gemeten door patentaanvragen. Een patentaanvraag geeft aan dat een bedrijf een stap neemt in de richting van de commercialisatie van het gepaten-teerde product, maar het is geen maat voor de daadwerkelijke introductie van het product op de markt, de innovatie. Om conclusies te trekken over de procycliciteit van innovaties gemeten door patentaanvragen zal verder onderzoek moeten worden gedaan naar de timing van patentaanvragen en daadwerkelijke productintroducties.

Dit onderzoek beschrijft en analyseert twee alternatieven op het Poissonmodel met QML-eigenschap zoals ge¨ımplementeerd in Fabrizio en Tsolmon (2014). Fabrizio en Tsolmon (2014) schonen hun data echter op een andere manier op dan in dit onderzoek is gedaan, waardoor het lastig is deze twee onderzoeken volledig met elkaar te vergelijken. Vervolgonderzoek zou zich kunnen richten op het schatten van het negatief binomiale en het zero-inflated Poisson-model met exact dezelfde dataset als gebruikt in Fabrizio en Tsolmon (2014). Een laatste suggestie voor vervolgonderzoek is analyseren of het zero-inflated negatief binomiale model (ZINB), zoals beschreven in Winkelmann (2008), beter bij de data past dan de besproken modellen in dit onderzoek. Het ZINB-model heeft mogelijk een variantiefunctie die beter bij de grote waarnemingen aansluit dan het zero-inflated Poissonmodel, waarbij variantie en ver-wachting van het gedeelte zonder additieve nullen gelijk worden verondersteld. Echter biedt de niet-significant van nul verschillende kans op additionele nullen uit het ZIP-model weinig perspectief voor het gebruik van een ander zero-inflated model.

7

Bibliografie

Aghion, P., Angeletos, G. M., Banerjee, A., Manova, K. (2010). Volatility and growth: Cre-dit constraints and the composition of investment. Journal of Monetary Economics, 57 (3), 246-265.

Aghion, P., Askenazy, P., Berman, N., Cette, G., Eymard, L. (2012). Credit Constraints and the Cyclicality of RD Investment: Evidence from France. Journal of the European Economic Association, 10 (5), 1001-1024.

Barlevy, G. (2007). On the Cyclicality of Research and Development. The American Econo-mic Review, 97 (4) 1131-1164.

Cohen, W. M., Nelson, R. R., Walsh, J. P. (2000). Protecting their Intellectual Assets: Appropriability Conditions and Why US Manufacturing Firms Patent (or not) (No. w7552). National Bureau of Economic Research.

Fabrizio, K. R., Tsolmon, U. (2014). An Empirical Examination of the Procyclicality of R&D Investment and Innovation. Review of Economics and Statistics, 96 (4), 662-675. Francois, P., Lloyd-Ellis, H. (2009). Schumpeterian Cycles with pro-Cyclical R&D. Re-view of Economic Dynamics, 12 (4), 567-591.

Griliches, Z. (1990). Patent Statistics as Economic Indicators: a Survey. National Bureau of Economic Research, 28, 1661-1707.

Heij, C., De Boer, P., Franses, P. H., Kloek, T., Van Dijk, H. K. (2004). Econometric Methods with Applications in Business and Economics. OUP Oxford.

Lambert, D. (1992). Zero-inflated Poisson Regression, with an Application to Defects in Manufacturing. Technometrics, 34 (1), 1-14.

McCullagh, P., Nelder, J.A. (1989). Generalized linear models, second edition. London: Chapman and Hall.

Mullahy, J. (1986). Specification and Testing of some Modified Count Data Models. Journal of Econometrics, 33 (3), 341-365.

Nelder, J.A., Wedderburn, R.W.M. (1972). Generalized linear models. Journal of Royal Statistical Society, series A, 135, 370-384.

of the Capitalist Process. New York: McGraw-Hill Book Company.

Winkelmann, R. (2008). Econometric Analysis of Count Data. Springer Science & Busi-ness Media.

8

Appendix

Appendix A

Afleiding van de gradi¨ent en Hessiaan van de log-likelihoodfunctie van het Poissonmodel. ∂L ∂λ = 1 + y λ ∂2L ∂λ2 = −y λ2 (8.0.1) Appendix B

Afleiding van de gradi¨ent en Hessiaan van de log-likelihoodfunctie van het negatief binomiale model. ∂L ∂λ = α α + y α α (α + λ)2 + y α + λ λ α + λ − λ (α + λ)2 = α(y − λ) λ(α + λ) ∂2L ∂λ2 = −αλ(α + λ) − α(y − λ)(α + 2λ) λ2_{(α + λ)}2 = αλ2− αy(α + 2λ) λ2_{(α + λ)}2

Voor het Fischerscorealgoritme wordt minus de verwachte waarde van de Hessiaan gebruikt als variantie-covariantiematrx (McCullagh & Nelder, 1989):

E ∂ L ∂λ = 0 en E ∂2L ∂λ2 = −α λ(α + λ)

Voor de schatter van de dispersieparameter, α, worden de eerste en tweede orde afgeleides

bepaald zoals hieronder, met DIG() en TRI() respectievelijk de digamma en trigamma functies, d.w.z. de eerste en tweede orde afgeleide van log(Γ()):

∂L ∂α = DIG(α + y) − DIG(α) + α α + λ α α + λ − α (α + λ)2 + log( α α + λ) + y α + λ λ −λ (α + λ)2

= DIG(α + y) − DIG(α) + log( α α + λ) + λ − y α + λ ∂2L ∂α2 = T RI(α + y) − T RI(α) + λ α(α + λ) − λ − y (α + λ)2 (8.0.2)

Appendix C

Afleiding van de gradi¨ent en Hessiaan van de log-likelihoodfunctie van het zero-inflated Pois-sonmodel.

Log-likelihoodfunctie: L(λ, ω) = (1 − I) ∗ log(ω + (1 − ω)e−λ) + I ∗ (log(1 − ω) − λ + ylog(λ))

Om ervoor te zorgen dat ω, de kans op y=0, tijdens de iteratieve gewogen kleinste kwadraten methode niet groter dan 1 wordt, wordt de getransformeerde variabele geschat: ξ = log( ω

1 − ω). ∂L ∂λ =    −(1−ω)e−λ ω+(1−ω)e−λ = −1 ω 1−ωeλ+1 = _1+e−1ξ+λ als y = 0 −1 + y_λ als y > 0 ∂2L ∂λ2 = ( _eλ+ξ (eλ+ξ₊₁₎2 = z(1 − z) als y = 0 met z = 1 1+eξ+λ −1 (1−ω)2 als y > 0 ∂L ∂ω = ( _1−e−λ ω+(1−ω)e−λ = e −1 ωeλ_(1−ω) = eλ₋₁ 1−ωz −1 1−ω ∂2L ∂ω2 = ( −(eλ_{− 1)ωe}λ_{+ (1 − ω)}−2 (eλ_{− 1) =} (eλ−1)2 (1−ω)2 z2 −1 (1−ω)2 ∂L ∂ξ = ∂L ∂ω ∂ω ∂ξ = ( eλ₋₁ 1−ωzω(1 − ω) = ω(e λ_{− 1)z} −1 1−ωω(1 − ω) = −ω ∂2L ∂ξ2 = ∂L ∂2_ω2 ∂2ω ∂ξ2 = (_(eλ₋₁₎2 (1−ω)2z2ω(1 − ω)(2(1 − ω) − 1) = (eλ₋₁₎2 1−ω z2ω(2(1 − ω) − 1) −1 (1−ω)2ω(1 − ω)(2(1 − ω) − 1) = −1 1−ωω(2(1 − ω) − 1)

De volledige gradi¨ent- en hessiaanfuncties worden vervolgens opgesteld met de volgende index functie: I =

(

0 als y = 0 1 als y > 0

Appendix D

R-script voor iteratief gewogen kleinste kwadraten voor het Poissonmodel met QML-eigenschap.

library("lfe")

family <- "quasipoisson"

# Read data and define year / gvkey as factor data <- cleandata

data$year <- as.factor(data$year) data$gvkey <- as.factor(data$gvkey) # Define covars; create modelformula

listCovars <- c("ln_lagrd" ,"ln_lag2rd" ,"ln_lagsales" ,"ln_lagemp" , "ln_lagppe" ,"ln_rgout" ,"ln_rgoutXpateff" ,

"ln_rgoutXfastcont" ,"ln_rgoutXext") model <- paste0(listCovars, collapse=" + ")

modelNB <- as.formula(paste0("NumPats ~ gvkey + year + ", model)) modelFELM <- as.formula(paste0("yadj ~ ", model, "+ year | gvkey")) # IRLS using felm()

# Using twice the log-likelihood as stopping criterion NumPats <- data$NumPats

devold <- 0

eta <- log(NumPats + 0.5*(NumPats==0)) mu <- exp(eta)

dev0 <- 2*sum(-NumPats + NumPats*log(NumPats + (NumPats==0))) for (iter in 1:50) {

# Fisher scoring for linear parameters data$weight <- mu

data$yadj <- eta + (NumPats-mu)/mu

felm <- felm(modelFELM, weight=data$weight, data=data) eta <- data$yadj - felm$residuals

mu <- exp(eta) # Check convergence

devFELM <- dev0 - 2*sum(-mu + NumPats*log(mu + (mu==0))) tol <- (abs(devFELM-devold)/devFELM)

if (TRUE) {

# Monitoring information if (iter == 1) {

cat("\n cycle deviance criterion\n") }

cat(paste0(formatC(iter, format="d", width=7),

formatC(devFELM, format="e", width=14, digits=4), formatC(tol, format="e", width=13, digits=2), "\n")) } if (tol < 1.0e-12) { break } devold <- devFELM } dfFELM <- felm$df

pearFELM <- sum( (NumPats-mu)*(NumPats-mu)/mu ) muPois <- mu

# Save regression output for qpoisson felmP <- felm

# Extract estimates and correct se estiFELM.P <- c(coef(felmP))

seFELM.P <- c(coef(summary(felmP))[, 2]/summary(felmP)$sigma) if (family == "quasipoisson") {

# glm() uses Pearson’s statistic to estimate the dispersion dispersion <- pearFELM/felmP$df

# seFELM.P <- sqrt(dispersion)*seFELM.P }

tvalFELM.P <- estiFELM.P/seFELM.P

pvalFELM.P <- round(10000*2*pnorm(-abs(tvalFELM.P)))/10000 cbind(estiFELM.P, seFELM.P, tvalFELM.P, pvalFELM.P)

# Pearson-test statistic

Appendix E

Fabrizio en Tsolmon (2014) gebruiken voor het analyseren van relatie (3.0.1) met het Pois-sonmodel met QML-eigenschap het pakket xtpqml in Stata, dat vaste effecten schat voor bedrijven. Dit pakket geeft de normale standaardafwijkingen uit het Poissonmodel, dus zon-der dat vermenigvuldigd wordt met de wortel van de dispersieparameter zoals in (2.2.4). In dit onderdeel wordt aangetoond dat het schatten van vergelijking (3.0.1) dezelfde resultaten oplevert met xtpqml in Stata als met de iteratief gewogen kleinste kwadratenmethode voor het Poissonmodel in R, wanneer dezelfde dataset wordt gebruikt. In de implementatie in R worden de standaardafwijkingen in dit geval ook niet vermenigdvuldigd met de wortel van de dispersieparameter, omdat dit in Stata ook niet gebeurd. Zoals hieronder te zien is, zijn de co¨effici¨entwaardes en standaardafwijkingen voor beide implementaties gelijk.

Implementatie in Stata

. *Fixed effect quasi poisson regression*

. xtpqml numpat ln_lagrd ln_lag2rd ln_lagsales ln_lagemp ln_lagppe ln_rgout ln_rgoutxpateff ln_rgoutxfastcont ln_rgoutxext yy*

> * , fe i(gvkey)

note: 127 groups (127 obs) dropped because of only one obs per group note: 748 groups (5931 obs) dropped because of all zero outcomes note: yy7 dropped because of collinearity

Iteration 0: log likelihood = -216518.57 Iteration 1: log likelihood = -142206.96 Iteration 2: log likelihood = -137632.07 Iteration 3: log likelihood = -137561.37 Iteration 4: log likelihood = -137561.28 Iteration 5: log likelihood = -137561.28

Conditional fixed-effects Poisson regression Number of obs = 42916

Group variable: gvkey Number of groups = 3569

Obs per group: min = 2

avg = 12.0

max = 28

Wald chi2(36) = 126945.99

Log likelihood = -137561.28 Prob > chi2 = 0.0000

---numpats | Coef. Std. Err. z P>|z| [95% Conf. Interval] ---+---ln_lagrd | .0285912 .0019942 14.34 0.000 .0246827 .0324997 ln_lag2rd | .021333 .0015642 13.64 0.000 .0182672 .0243988 ln_lagsales | .1849698 .0044375 41.68 0.000 .1762724 .1936672 ln_lagemp | .1961478 .0038238 51.30 0.000 .1886533 .2036424 ln_lagppe | .0879749 .0044721 19.67 0.000 .0792098 .09674 ln_rgout | .835678 .0186114 44.90 0.000 .7992003 .8721556 ln_rgoutxpateff | -1.481747 .0384443 -38.54 0.000 -1.557096 -1.406398 ln_rgoutxfastcont | -.0188446 .0077656 -2.43 0.015 -.0340648 -.0036244 ln_rgoutxext | .0943006 .0042191 22.35 0.000 .0860313 .1025698 yy1 | .2708326 .0117003 23.15 0.000 .2479004 .2937648 yy2 | .2575761 .0116227 22.16 0.000 .234796 .2803561 yy3 | .2126692 .0115832 18.36 0.000 .1899666 .2353717 yy4 | .1224005 .0116329 10.52 0.000 .0996005 .1452005 yy5 | .0678607 .0116822 5.81 0.000 .044964 .0907574 yy6 | .0591932 .0115696 5.12 0.000 .0365173 .0818691 yy8 | .0321368 .0115877 2.77 0.006 .0094254 .0548482 yy9 | -.0167139 .0117305 -1.42 0.154 -.0397053 .0062775 yy10 | .0392519 .0114742 3.42 0.001 .0167628 .061741 yy11 | .0703978 .0113454 6.20 0.000 .0481612 .0926343 yy12 | .0976867 .0113237 8.63 0.000 .0754926 .1198808 yy13 | .1211658 .0112305 10.79 0.000 .0991544 .1431772 yy14 | .1735549 .0110151 15.76 0.000 .1519657 .1951442 yy15 | .2127118 .0108654 19.58 0.000 .1914159 .2340076 yy16 | .2111239 .0108266 19.50 0.000 .1899041 .2323436 yy17 | .1804151 .0108011 16.70 0.000 .1592454 .2015847 yy18 | .2024514 .0108173 18.72 0.000 .1812499 .2236529 yy19 | .2103888 .010759 19.55 0.000 .1893015 .2314761 yy20 | .2943892 .0105964 27.78 0.000 .2736207 .3151577 yy21 | .4195228 .0103362 40.59 0.000 .3992641 .4397814 yy22 | .3030874 .0104785 28.92 0.000 .2825498 .3236249 yy23 | .4220758 .0103842 40.65 0.000 .4017233 .4424284 yy24 | .337229 .0105208 32.05 0.000 .3166086 .3578493 yy25 | .3384032 .0106399 31.81 0.000 .3175494 .359257 yy26 | .3328651 .0107231 31.04 0.000 .3118482 .353882 yy27 | .2903733 .0107849 26.92 0.000 .2692352 .3115114 yy28 | .2081156 .0109395 19.02 0.000 .1866747 .2295565

> library("lfe") >

> family <- "quasipoisson" >

> # Read data and define year / gvkey as factor > data <- cleandata

> data$year <- as.factor(data$year) > data$gvkey <- as.factor(data$gvkey) >

> # Define covars; create modelformula

> listCovars <- c("ln_lagrd" ,"ln_lag2rd" ,"ln_lagsales" ,"ln_lagemp" , + "ln_lagppe" ,"ln_rgout" ,"ln_rgoutXpateff" ,

+ "ln_rgoutXfastcont" ,"ln_rgoutXext") > model <- paste0(listCovars, collapse=" + ")

> modelNB <- as.formula(paste0("NumPats ~ gvkey + year + ", model)) > modelFELM <- as.formula(paste0("yadj ~ ", model, "+ year | gvkey")) >

> # IRLS using felm()

> # Using twice the log-likelihood as stopping criterion > NumPats <- data$NumPats

> devold <- 0

> eta <- log(NumPats + 0.5*(NumPats==0)) > mu <- exp(eta)

> dev0 <- 2*sum(-NumPats + NumPats*log(NumPats + (NumPats==0))) > for (iter in 1:50) {

+ # Fisher scoring for linear parameters + data$weight <- mu

+ data$yadj <- eta + (NumPats-mu)/mu

+ felm <- felm(modelFELM, weight=data$weight, data=data) + eta <- data$yadj - felm$residuals

+ mu <- exp(eta) + # Check convergence

+ devFELM <- dev0 - 2*sum(-mu + NumPats*log(mu + (mu==0))) + tol <- (abs(devFELM-devold)/devFELM)

+ if (TRUE) {

+ # Monitoring information + if (iter == 1) {

+ cat("\n cycle deviance criterion\n")

+ }

+ cat(paste0(formatC(iter, format="d", width=7),

+ formatC(tol, format="e", width=13, digits=2), "\n")) + } + if (tol < 1.0e-12) { + break + } + devold <- devFELM + }

cycle deviance criterion 1 2.4539e+05 1.00e+00 2 2.0595e+05 1.92e-01 3 2.0258e+05 1.66e-02 4 2.0227e+05 1.51e-03 5 2.0220e+05 3.71e-04 6 2.0217e+05 1.35e-04 7 2.0216e+05 4.96e-05 8 2.0216e+05 1.82e-05 9 2.0216e+05 6.71e-06 10 2.0216e+05 2.47e-06 11 2.0216e+05 9.08e-07 12 2.0216e+05 3.34e-07 13 2.0216e+05 1.23e-07 14 2.0216e+05 4.52e-08 15 2.0216e+05 1.66e-08 16 2.0216e+05 6.12e-09 17 2.0216e+05 2.25e-09 18 2.0216e+05 8.28e-10 19 2.0216e+05 3.05e-10 20 2.0216e+05 1.12e-10 21 2.0216e+05 4.12e-11 22 2.0216e+05 1.52e-11 23 2.0216e+05 5.58e-12 24 2.0216e+05 2.05e-12 25 2.0216e+05 7.56e-13 > dfFELM <- felm$df

> pearFELM <- sum( (NumPats-mu)*(NumPats-mu)/mu ) >

> # Save regression output for qpoisson > felmP <- felm

>

> estiFELM.P <- c(coef(felmP))

> seFELM.P <- c(coef(summary(felmP))[, 2]/summary(felmP)$sigma) > if (family == "quasipoisson") {

+ # glm() uses Pearson’s statistic to estimate the dispersion + dispersion <- pearFELM/felmP$df

+ # seFELM.P <- sqrt(dispersion)*seFELM.P + }

> tvalFELM.P <- estiFELM.P/seFELM.P

> pvalFELM.P <- round(10000*2*pnorm(-abs(tvalFELM.P)))/10000 > cbind(estiFELM.P, seFELM.P, tvalFELM.P, pvalFELM.P)

estiFELM.P seFELM.P tvalFELM.P pvalFELM.P ln_lagrd 0.02859121 0.001994186 14.337284 0.0000 ln_lag2rd 0.02133300 0.001564216 13.638144 0.0000 ln_lagsales 0.18496978 0.004437542 41.682935 0.0000 ln_lagemp 0.19614780 0.003823822 51.296271 0.0000 ln_lagppe 0.08797490 0.004472095 19.671968 0.0000 ln_rgout 0.83567797 0.018611397 44.901411 0.0000 ln_rgoutXpateff -1.48174708 0.038444262 -38.542737 0.0000 ln_rgoutXfastcont -0.01884462 0.007765565 -2.426690 0.0152 ln_rgoutXext 0.09430059 0.004219082 22.350974 0.0000 year1976 -0.01325654 0.011342772 -1.168722 0.2425 year1977 -0.05816345 0.011428574 -5.089301 0.0000 year1978 -0.14843213 0.011599014 -12.796961 0.0000 year1979 -0.20297191 0.011726716 -17.308503 0.0000 year1980 -0.21163942 0.011577195 -18.280716 0.0000 year1981 -0.27083262 0.011700326 -23.147443 0.0000 year1982 -0.23869585 0.011552930 -20.661065 0.0000 year1983 -0.28754652 0.011719383 -24.535977 0.0000 year1984 -0.23158071 0.011531115 -20.083115 0.0000 year1985 -0.20043487 0.011412890 -17.562149 0.0000 year1986 -0.17314592 0.011402683 -15.184666 0.0000 year1987 -0.14966684 0.011370359 -13.162895 0.0000 year1988 -0.09727771 0.011211593 -8.676529 0.0000 year1989 -0.05812086 0.011082527 -5.244369 0.0000 year1990 -0.05970876 0.011039980 -5.408412 0.0000 year1991 -0.09041756 0.011002825 -8.217668 0.0000 year1992 -0.06838121 0.011058107 -6.183808 0.0000 year1993 -0.06044381 0.011022550 -5.483650 0.0000 year1994 0.02355658 0.010898078 2.161535 0.0307 year1995 0.14869015 0.010681881 13.919847 0.0000 year1996 0.03225474 0.010857297 2.970789 0.0030

year1997 0.15124322 0.010814258 13.985539 0.0000 year1998 0.06639637 0.010974603 6.050002 0.0000 year1999 0.06757055 0.011117351 6.077936 0.0000 year2000 0.06203247 0.011235560 5.521084 0.0000 year2001 0.01954070 0.011315503 1.726896 0.0842 year2002 -0.06271702 0.011471491 -5.467207 0.0000

Appendix F

De R-functie glm kan een maximumlikelihoodschatting uitvoeren voor een Poissonmodel met een kleine hoeveelheid vaste effecten. Wanneer de hoeveelheid vaste effecten te groot wordt, zoals het geval is in dit onderzoek met 4444 bedrijven als vaste effecten, convergeert deze functie niet en moet de iteratief gewogen kleinste kwadratenmethode gebruikt worden. In dit onderdeel wordt aangetoond dat het schatten van vergelijking (3.0.1) met het Poissonmo-del met QML-eigenschap dezelfde resultaten oplevert voor de R-functie glm als voor itera-tief gewogen kleinste kwadraten. In het onderstaande uitgevoerde R-script is te zien dat de co¨effici¨entwaardes, standaardafwijkingen voor beide implementaties gelijk zijn. Onderaan is te zien dat de deviance (devGLM en devFELM), dispersieparameter (pearGLM en pearFELM) en het aantal vrijheidsgraden (dfGLM en dfFELM) hetzelfde zijn voor beide implementaties. > family <- "quasipoisson"

> library("lfe") >

> # Read data and define year / gvkey as factor > data <- cleandata

> data$year <- as.factor(data$year) > data$gvkey <- as.factor(data$gvkey) >

> if (TRUE) {

+ # Select 10 random companies + nsites <- 10

+ set.seed(873021873)

+ levels <- levels(data$gvkey) + length(levels)

+ subsample <- sample(levels, nsites)

+ data <- data[(data$gvkey %in% subsample),] + data$year <- as.factor(as.numeric(data$year)) + data$gvkey <- as.factor(as.numeric(data$gvkey)) + }

>

> # Define covars; create modelformulae

> # Note that gvkey must be in front for modelNB/modelLM

> listCovars <- c("ln_lagrd" ,"ln_lag2rd" ,"ln_lagsales" ,"ln_lagemp" , + "ln_lagppe" ,"ln_rgout" ,"ln_rgoutXpateff" ,

+ "ln_rgoutXfastcont" ,"ln_rgoutXext") > model <- paste0(listCovars, collapse=" + ")

> modelNB <- as.formula(paste0("NumPats ~ ", model, " + gvkey + year")) > modelFELM <- as.formula(paste0("yadj ~ ", model, " + year | gvkey"))

> # Fit Poisson regression; save parameter estimates for comparison > glm <- glm(modelNB, data=data, family=family)

> summary(glm) Call:

glm(formula = modelNB, family = family, data = data) Deviance Residuals:

Min 1Q Median 3Q Max

-3.7172 -1.0394 -0.2597 0.5626 4.0667 Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) -2.33569 3.74411 -0.624 0.534091 ln_lagrd 0.37427 0.19475 1.922 0.057339 . ln_lag2rd -0.07199 0.18289 -0.394 0.694668 ln_lagsales -0.38882 0.32213 -1.207 0.230136 ln_lagemp -0.03189 0.31444 -0.101 0.919411 ln_lagppe 0.88300 0.23250 3.798 0.000245 *** ln_rgout 2.82559 7.52317 0.376 0.707983 ln_rgoutXpateff -15.15523 12.64881 -1.198 0.233554 ln_rgoutXfastcont 1.27748 1.21136 1.055 0.294037 ln_rgoutXext -0.85482 0.97310 -0.878 0.381704 gvkey161 17.07404 9.47675 1.802 0.074466 . gvkey1065 -14.46125 1653.32424 -0.009 0.993038 gvkey1244 3.38106 0.79914 4.231 4.99e-05 *** gvkey1952 -9.08195 21.35217 -0.425 0.671459 gvkey1995 -0.11401 2.70369 -0.042 0.966445 gvkey2333 51.30065 25.93452 1.978 0.050539 . gvkey2862 -2.09157 3.50493 -0.597 0.551957 gvkey4058 -1.07469 4.71753 -0.228 0.820240 gvkey4220 -1.65350 4.30872 -0.384 0.701935 year2 -0.30249 0.18826 -1.607 0.111103 year3 -0.35468 0.18420 -1.926 0.056866 . year4 -0.31579 0.19350 -1.632 0.105675 year5 -0.09471 0.20567 -0.460 0.646112 year6 -0.46487 0.23074 -2.015 0.046491 * year7 -0.33482 0.22338 -1.499 0.136901 year8 -0.46150 0.23103 -1.998 0.048345 * year9 -0.53100 0.24130 -2.201 0.029957 * year10 -0.54461 0.25376 -2.146 0.034162 *

year11 -0.65387 0.26634 -2.455 0.015729 * year12 -0.39473 0.26918 -1.466 0.145533 year13 -0.37612 0.28135 -1.337 0.184160 year14 -0.32464 0.27304 -1.189 0.237127 year15 -0.36613 0.28129 -1.302 0.195909 year16 -0.33317 0.29230 -1.140 0.256952 year17 -0.24877 0.30037 -0.828 0.409422 year18 -0.36120 0.30677 -1.177 0.241696 year19 -0.31974 0.28883 -1.107 0.270814 year20 -0.27777 0.27693 -1.003 0.318152 year21 -0.33559 0.26140 -1.284 0.202038 year22 -0.28749 0.25166 -1.142 0.255906 year23 -0.29141 0.26892 -1.084 0.281003 year24 -0.40937 0.29670 -1.380 0.170596 year25 -0.58575 0.30118 -1.945 0.054468 . year26 -0.77655 0.31922 -2.433 0.016680 * year27 -0.76486 0.31844 -2.402 0.018067 * year28 -0.83891 0.30084 -2.789 0.006287 **

---Signif. codes: 0 ^a***^a 0.001 ^a**^a 0.01 ^a*^a 0.05 ^a.^a 0.1 ^a ^a 1 (Dispersion parameter for quasipoisson family taken to be 3.703141)

Null deviance: 21367.03 on 150 degrees of freedom Residual deviance: 293.86 on 105 degrees of freedom AIC: NA

Number of Fisher Scoring iterations: 15 > estiGLM <- coef(glm)

> seGLM <- coef(summary(glm))[, 2] > thetaGLM <- glm$theta

> thetaseGLM <- glm$SE.theta > devGLM <- glm$deviance

> pearGLM <- sum(resid(glm,type="pear") * resid(glm,type="pear")) > dfGLM <- glm$df.residual

>

> # IRLS using felm()

> # Using twice the log-likelihood as stopping criterion > NumPats <- data$NumPats

> eta <- log(NumPats + 0.5*(NumPats==0)) > mu <- exp(eta)

> dev0 <- 2*sum(-NumPats + NumPats*log(NumPats + (NumPats==0))) > for (iter in 1:50) {

+ # Fisher scoring for linear parameters + data$weight <- mu

+ data$yadj <- eta + (NumPats-mu)/mu

+ felm <- felm(modelFELM, weight=data$weight, data=data) + eta <- data$yadj - felm$residuals

+ mu <- exp(eta) + # Check convergence

+ devFELM <- dev0 - 2*sum(-mu + NumPats*log(mu + (mu==0))) + tol <- (abs(devFELM-devold)/devFELM)

+ if (TRUE) {

+ # Monitoring information + if (iter == 1) {

+ cat("\n cycle deviance criterion\n")

+ }

+ cat(paste0(formatC(iter, format="d", width=7),

+ formatC(devFELM, format="e", width=14, digits=4), + formatC(tol, format="e", width=13, digits=2), "\n"))

+ } + if (tol < 1.0e-12) { + break + } + devold <- devFELM + }

cycle deviance criterion 1 3.3896e+02 1.00e+00 2 2.9773e+02 1.39e-01 3 2.9464e+02 1.05e-02 4 2.9414e+02 1.73e-03 5 2.9396e+02 5.99e-04 6 2.9389e+02 2.20e-04 7 2.9387e+02 8.11e-05 8 2.9386e+02 2.98e-05 9 2.9386e+02 1.10e-05 10 2.9386e+02 4.04e-06 11 2.9386e+02 1.49e-06 12 2.9386e+02 5.47e-07

13 2.9386e+02 2.01e-07 14 2.9386e+02 7.40e-08 15 2.9386e+02 2.72e-08 16 2.9386e+02 1.00e-08 17 2.9386e+02 3.68e-09 18 2.9386e+02 1.35e-09 19 2.9386e+02 4.98e-10 20 2.9386e+02 1.83e-10 21 2.9386e+02 6.74e-11 22 2.9386e+02 2.48e-11 23 2.9386e+02 9.16e-12 24 2.9386e+02 3.32e-12 25 2.9386e+02 1.24e-12 26 2.9386e+02 4.46e-13 > dfFELM <- felm$df

> pearFELM <- sum( (NumPats-mu)*(NumPats-mu)/mu ) >

> # Extract estimates and correct se > estiFELM <- c(coef(felm))

> seFELM <- c(coef(summary(felm))[, 2]/summary(felm)$sigma) > if (family == "quasipoisson") {

+ # glm() uses Pearson’s statistic to estimate the dispersion + dispersion <- pearFELM/felm$df

+ seFELM <- sqrt(dispersion)*seFELM + }

> tvalFELM <- estiFELM/seFELM

> pvalFELM <- round(10000*2*pnorm(-abs(tvalFELM)))/10000 > cbind(estiFELM, seFELM, tvalFELM, pvalFELM)

estiFELM seFELM tvalFELM pvalFELM ln_lagrd 0.37427309 0.1947493 1.9218202 0.0546 ln_lag2rd -0.07198818 0.1828925 -0.3936092 0.6939 ln_lagsales -0.38881581 0.3221280 -1.2070226 0.2274 ln_lagemp -0.03189011 0.3144389 -0.1014191 0.9192 ln_lagppe 0.88300193 0.2324991 3.7978730 0.0001 ln_rgout 2.82559134 7.5231695 0.3755852 0.7072 ln_rgoutXpateff -15.15523309 12.6488053 -1.1981553 0.2309 ln_rgoutXfastcont 1.27748307 1.2113628 1.0545834 0.2916 ln_rgoutXext -0.85481762 0.9730958 -0.8784517 0.3797 year2 -0.30249278 0.1882586 -1.6067937 0.1081 year3 -0.35468114 0.1841978 -1.9255452 0.0542 year4 -0.31578624 0.1934961 -1.6320028 0.1027

year5 -0.09471249 0.2056748 -0.4604964 0.6452 year6 -0.46486758 0.2307364 -2.0147125 0.0439 year7 -0.33481635 0.2233753 -1.4988963 0.1339 year8 -0.46150031 0.2310262 -1.9976104 0.0458 year9 -0.53100323 0.2413026 -2.2005695 0.0278 year10 -0.54461409 0.2537597 -2.1461801 0.0319 year11 -0.65386892 0.2663363 -2.4550501 0.0141 year12 -0.39472589 0.2691839 -1.4663800 0.1425 year13 -0.37612081 0.2813477 -1.3368540 0.1813 year14 -0.32463865 0.2730381 -1.1889867 0.2344 year15 -0.36612537 0.2812934 -1.3015782 0.1931 year16 -0.33317488 0.2923031 -1.1398268 0.2544 year17 -0.24877391 0.3003700 -0.8282249 0.4075 year18 -0.36119578 0.3067731 -1.1774036 0.2390 year19 -0.31973596 0.2888254 -1.1070214 0.2683 year20 -0.27777105 0.2769318 -1.0030305 0.3158 year21 -0.33559077 0.2614041 -1.2838005 0.1992 year22 -0.28748606 0.2516614 -1.1423525 0.2533 year23 -0.29140968 0.2689153 -1.0836484 0.2785 year24 -0.40936718 0.2966976 -1.3797455 0.1677 year25 -0.58575055 0.3011803 -1.9448500 0.0518 year26 -0.77655085 0.3192216 -2.4326386 0.0150 year27 -0.76485695 0.3184373 -2.4019079 0.0163 year28 -0.83891348 0.3008436 -2.7885367 0.0053 >

> # Compare glm() and iterative weightes least squares > cbind(devGLM, devFELM, pearGLM, pearFELM, dfGLM, dfFELM)

devGLM devFELM pearGLM pearFELM dfGLM dfFELM [1,] 293.8564 293.8564 388.8298 388.8298 105 105