In [40]:
import pandas as pd
from pandas.tools.plotting import lag_plot
from pandas.tools.plotting import autocorrelation_plot
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from matplotlib.ticker import FuncFormatter
import math
import os
import time
path = '.'
saveFigs = True
In [41]:
# Find the data files, takes the most recent if multiple
tests = ['test1a', 'test1b', 'test2a', 'test2b', 'test3a', 'test3b', 'test4a', 'test4b']
testdatafiles = []
fail = False
for test in tests:
try:
testdatafiles.append(sorted((fn for fn in os.listdir(path) if test in fn and 'rundata' in fn))[0])
except:
print "no data for "+test
fail = True
if fail==False:
test1A = pd.DataFrame.from_csv(testdatafiles[0])
test1B = pd.DataFrame.from_csv(testdatafiles[1])
test2A = pd.DataFrame.from_csv(testdatafiles[2])
test2B = pd.DataFrame.from_csv(testdatafiles[3])
test3A = pd.DataFrame.from_csv(testdatafiles[4])
test3B = pd.DataFrame.from_csv(testdatafiles[5])
test4A = pd.DataFrame.from_csv(testdatafiles[6])
test4B = pd.DataFrame.from_csv(testdatafiles[7])
typdatafiles = sorted((fn for fn in os.listdir(path) if 'typical' in fn))
typRunFile = next(x for x in typdatafiles if 'run' in x)
typRun = pd.DataFrame.from_csv(path+'/'+typRunFile)
typInnovFile = next(x for x in typdatafiles if 'innov' in x)
typInnov = pd.DataFrame.from_csv(path+'/'+typInnovFile)
typStepFile = next(x for x in typdatafiles if 'step' in x)
typStep = pd.DataFrame.from_csv(path+'/'+typStepFile)
fileTime = time.strftime("%Y-%m-%d %H%M_%S")
testPatentLife = [5, 200]
testMonopProfit = [1, 3, 5]
radCutOff = 2
numSteps = len(typStep[typStep.run==0])
In [42]:
# Generic plotting method for maximum BPF reached
def doMaxBPFPlot(y1, y2, y3, x, title1, title2, xlab, invertx):
# Define fitted lines
zn = np.polyfit(x, y1, 4)
pn = np.poly1d(zn)
z5 = np.polyfit(x,y2,4)
p5 = np.poly1d(z5)
z200 = np.polyfit(x,y3,4)
p200 = np.poly1d(z200)
# Create x axis variable
x2 = np.linspace(round(x.values[0],3), round(x.values[-1],3), 10)
# Create plots
fig = plt.figure(figsize=(13,5))
ax1 = fig.add_subplot(121)
ax1.plot(x2, pn(x2), lw=2)
ax1.plot(x2, p5(x2), lw=2, ls='--', c='r')
ax1.plot(x2, p200(x2), lw=2, ls=':', c='cyan')
if invertx:
ax1.invert_xaxis()
ax1.set_title(title1)
ax1.set_ylabel('Max BPF Height')
ax1.set_xlabel(xlab)
labels = ['No Patents', 'Low Duration', 'High Duration']
ax1.legend(labels, loc=0)
ax2 = fig.add_subplot(122)
ax2.plot(x2, (p5(x2)-pn(x2))/((pn(x2)+p5(x2))/2), lw=2, ls='--', c='r')
ax2.plot(x2, (p200(x2)-pn(x2))/((pn(x2)+p200(x2))/2), lw=2, ls=':', c='cyan')
ax2.axhline(0, color='black')
if invertx:
ax2.invert_xaxis()
ax2.set_title(title2)
ax2.set_ylabel('Percent Diff (Pat-NoPat)')
ax2.set_xlabel(xlab)
labels = ['Low Duration', 'High Duration']
ax2.legend(labels, loc=0)
return fig
In [43]:
# Generate example profit decay figure
monopProfit = 1
piDecayRate = .1
profitStream0 = []
for i in range(10):
profitStream0.append(round(monopProfit * math.exp(-.1*i),3))
profitStream1 = []
for i in range(10):
profitStream1.append(round(monopProfit * math.exp(-.2*i),3))
profitStream2 = []
for i in range(10):
profitStream2.append(round(monopProfit * math.exp(-.3*i),3))
print profitStream0
x = range(0,len(profitStream0))
fig = plt.figure(figsize=(5,5))
ax1 = fig.add_subplot(111)
ax1.plot(x, profitStream0, lw=2, c='k')
ax1.plot(profitStream1, lw=2, ls='--', c='b')
ax1.plot(profitStream2, lw=2, ls='-.', c='r')
labels = ['Profit Decay=0.1', 'Profit Decay=0.2', 'Profit Decay=0.3']
ax1.legend(labels)
ax1.set_title('Example Profit Decay Rates')
ax1.set_ylabel('Monopoly Rents')
ax1.set_xlabel('Time')
if saveFigs:
fileName = "profDecayExample.png"
plt.savefig(fileName)
plt.show()
The parameters for a typical model run, without patents, are shown in Appendix B (of paper). A histogram of innovation sizes, captured during the simulation, is shown in the figure below. The distribution of innovation sizes is highly skewed, agreeing with the empirical observations described above. Most innovations increased the BPF for the given column by only one row. A log-log plot of innovation sizes is shown in the figure below. The log-log plot shows a roughly linear decrease in frequency, suggesting a power-law or Poisson process. A simple ordinary least squares regression on the log-log plot gives a coefficient of -0.91.
In [44]:
# Select the innovation sizes for one of the runs
innovSizes = typInnov.innovSize
# Plot histogram of innovation size
fig = plt.figure(figsize=(13,5))
ax1 = fig.add_subplot(121)
ax1.hist(innovSizes, bins=20)
plt.xlabel('Innovation Size')
plt.ylabel('Frequency')
plt.title('Histogram of Innovation Size')
# Create bins for histogram
inputBins = [0, 1.5]
maxInnovSize = max(innovSizes)
while inputBins[-1] < maxInnovSize:
inputBins.append(inputBins[-1]*1.5)
# Generate histogram
hist, bins = np.histogram(innovSizes, bins=inputBins)
# Create log series for log plots, replacing negative infinity values with zero
bins = list(inputBins)
bins.remove(0)
lbins = list(np.log(bins))
lhist = list(np.log(hist))
for i in range(len(lhist)):
if lhist[i] == float('-inf'):
lhist[i] = 0
# Fit first degree polynomial to log series
coeff = np.polyfit(lbins,lhist,1)
polynomial = np.poly1d(coeff)
ys = polynomial(lbins)
logCoef = round(polynomial[1],2)
# Plot log series
ax2 = fig.add_subplot(122)
ax2.plot(lbins,lhist)
ax2.plot(lbins, ys)
plt.xlabel('Log Innovation Size')
plt.ylabel('Log Frequency')
plt.title('Log Plot of Innovation Size')
plt.text(max(lbins)*0.6, max(ys)*0.8, 'OLS Coefficient: {0}'.format(round(polynomial[1],2)))
if saveFigs:
figName = "{0}_TypicalRun_InnovHist.png".format(fileTime)
plt.savefig(figName)
plt.show()
A time series plot of the change in BPF height each step is shown in the figure below. This plot captures the sum of all innovation sizes for each step. Most steps see little change in BPF height while very few see massive increases. A lag-plot of the change in BPF height each step is also shown in the figure below. This lag-plot shows that the BPF data does not exhibit very much autocorrelation.
In [45]:
# Select BPF Change time series for one of the runs
BPFChange = typStep[typStep.run==0].BPFChange
# Plot BPF Change time series
fig = plt.figure(figsize=(13,5))
ax1 = fig.add_subplot(121)
ax1.plot(BPFChange)
plt.xlabel('Step')
plt.ylabel('BPF Change')
plt.title('BPF Change Time Series')
# BPF Change lag plot
ax1 = fig.add_subplot(122)
lag_plot(BPFChange)
plt.xlabel('BPFChange(t)')
plt.ylabel('BPFChange(t + 1)')
plt.title('BPF Change Lag Plot')
if saveFigs:
figName = "{0}_TypicalRun_BPFChange.png".format(fileTime)
plt.savefig(figName)
plt.show()
To test for temporal clustering in the arrival of radical innovations a time series of counts of radical innovations per step is generated and analyzed. Following the analysis of Silverberg and Verspagen 2003, Poisson and negative binomial regressions are used to determine if the data exhibit over-dispersion, i.e. variance greater than that expected by the Poisson model. Because the definition of radical with respect to the percolation model output is even more arbitrary, several definitions of radical are tested.
In [46]:
def pois_nbin_regs(endog, exog):
# Run Possion regression
poisson_mod = sm.Poisson(endog, exog)
poisson_res = poisson_mod.fit(method="newton")
poisson_sum = poisson_res.summary()
# Run Negative Binomial regression
nbin_mod = sm.NegativeBinomial(endog, exog)
nbin_res = nbin_mod.fit(disp=False)
nbin_sum = nbin_res.summary()
return poisson_sum, poisson_res, nbin_sum, nbin_res
def doRegs(rad, outputDF):
radCutOff=rad
# Create a list of radical innovation counts, one for each step
innovCounts = []
for i in range(numSteps):
innovSizeDF = typInnov[(typInnov.run==0) & (typInnov.step==i)].innovSize
if len(innovSizeDF)>0:
radCounter = 0
for j in innovSizeDF.index:
if innovSizeDF[j] >=radCutOff:
radCounter += 1
innovCounts.append(radCounter)
else:
innovCounts.append(0)
# Create endogenous and exogenous variables
endog = pd.DataFrame(innovCounts, index=range(0,len(innovCounts)), columns=['innovCounts'])
exog = pd.DataFrame(index=range(0,numSteps))
exog = sm.add_constant(exog, prepend = False)
poisson_sum, poisson_res, nbin_sum, nbin_res = pois_nbin_regs(endog,exog)
row = pd.Series(["Poisson0", radCutOff, poisson_sum.tables[0][0][1], poisson_sum.tables[0][0][3], poisson_sum.tables[0][1][1], poisson_sum.tables[0][1][3], poisson_sum.tables[0][2][1], poisson_sum.tables[0][2][3], poisson_sum.tables[0][3][1], poisson_sum.tables[0][3][3], poisson_sum.tables[0][4][1], poisson_sum.tables[0][4][3], poisson_sum.tables[0][5][1], poisson_sum.tables[0][5][3], poisson_sum.tables[0][6][3], poisson_res.params['const'],poisson_res.pvalues['const'],None,None,None,None,None,None,None,None],index=modelOutputCols)
outputDF = outputDF.append(row, ignore_index=True)
row = pd.Series(["NegBin0",radCutOff, nbin_sum.tables[0][0][1], nbin_sum.tables[0][0][3], nbin_sum.tables[0][1][1], nbin_sum.tables[0][1][3], nbin_sum.tables[0][2][1], nbin_sum.tables[0][2][3], nbin_sum.tables[0][3][1], nbin_sum.tables[0][3][3], nbin_sum.tables[0][4][1], nbin_sum.tables[0][4][3], nbin_sum.tables[0][5][1], nbin_sum.tables[0][5][3], nbin_sum.tables[0][6][3], nbin_res.params['const'],nbin_res.pvalues['const'],None,None,None,None,None,None,nbin_res.params['alpha'],nbin_res.pvalues['alpha']],index=modelOutputCols)
outputDF = outputDF.append(row, ignore_index=True)
exog['step'] = pd.Series(range(0,numSteps), index=exog.index)
poisson_sum, poisson_res, nbin_sum, nbin_res = pois_nbin_regs(endog,exog)
row = pd.Series(["Poisson1", radCutOff, poisson_sum.tables[0][0][1], poisson_sum.tables[0][0][3], poisson_sum.tables[0][1][1], poisson_sum.tables[0][1][3], poisson_sum.tables[0][2][1], poisson_sum.tables[0][2][3], poisson_sum.tables[0][3][1], poisson_sum.tables[0][3][3], poisson_sum.tables[0][4][1], poisson_sum.tables[0][4][3], poisson_sum.tables[0][5][1], poisson_sum.tables[0][5][3], poisson_sum.tables[0][6][3], poisson_res.params['const'],poisson_res.pvalues['const'],poisson_res.params['step'],poisson_res.pvalues['step'],None,None,None,None,None,None],index=modelOutputCols)
outputDF = outputDF.append(row, ignore_index=True)
row = pd.Series(["NegBin1",radCutOff, nbin_sum.tables[0][0][1], nbin_sum.tables[0][0][3], nbin_sum.tables[0][1][1], nbin_sum.tables[0][1][3], nbin_sum.tables[0][2][1], nbin_sum.tables[0][2][3], nbin_sum.tables[0][3][1], nbin_sum.tables[0][3][3], nbin_sum.tables[0][4][1], nbin_sum.tables[0][4][3], nbin_sum.tables[0][5][1], nbin_sum.tables[0][5][3], nbin_sum.tables[0][6][3], nbin_res.params['const'],nbin_res.pvalues['const'],nbin_res.params['step'],nbin_res.pvalues['step'],None,None,None,None,nbin_res.params['alpha'],nbin_res.pvalues['alpha']],index=modelOutputCols)
outputDF = outputDF.append(row, ignore_index=True)
exog['step2'] = pd.Series(exog.step**2, index=exog.index)
poisson_sum, poisson_res, nbin_sum, nbin_res = pois_nbin_regs(endog,exog)
row = pd.Series(["Poisson2", radCutOff, poisson_sum.tables[0][0][1], poisson_sum.tables[0][0][3], poisson_sum.tables[0][1][1], poisson_sum.tables[0][1][3], poisson_sum.tables[0][2][1], poisson_sum.tables[0][2][3], poisson_sum.tables[0][3][1], poisson_sum.tables[0][3][3], poisson_sum.tables[0][4][1], poisson_sum.tables[0][4][3], poisson_sum.tables[0][5][1], poisson_sum.tables[0][5][3], poisson_sum.tables[0][6][3], poisson_res.params['const'],poisson_res.pvalues['const'],poisson_res.params['step'],poisson_res.pvalues['step'],poisson_res.params['step2'],poisson_res.pvalues['step2'],None,None,None,None],index=modelOutputCols)
outputDF = outputDF.append(row, ignore_index=True)
row = pd.Series(["NegBin2",radCutOff, nbin_sum.tables[0][0][1], nbin_sum.tables[0][0][3], nbin_sum.tables[0][1][1], nbin_sum.tables[0][1][3], nbin_sum.tables[0][2][1], nbin_sum.tables[0][2][3], nbin_sum.tables[0][3][1], nbin_sum.tables[0][3][3], nbin_sum.tables[0][4][1], nbin_sum.tables[0][4][3], nbin_sum.tables[0][5][1], nbin_sum.tables[0][5][3], nbin_sum.tables[0][6][3], nbin_res.params['const'],nbin_res.pvalues['const'],nbin_res.params['step'],nbin_res.pvalues['step'],nbin_res.params['step2'],nbin_res.pvalues['step2'],None,None,nbin_res.params['alpha'],nbin_res.pvalues['alpha']],index=modelOutputCols)
outputDF = outputDF.append(row, ignore_index=True)
exog['step3'] = pd.Series(exog.step**3, index=exog.index)
poisson_sum, poisson_res, nbin_sum, nbin_res = pois_nbin_regs(endog,exog)
row = pd.Series(["Poisson3", radCutOff, poisson_sum.tables[0][0][1], poisson_sum.tables[0][0][3], poisson_sum.tables[0][1][1], poisson_sum.tables[0][1][3], poisson_sum.tables[0][2][1], poisson_sum.tables[0][2][3], poisson_sum.tables[0][3][1], poisson_sum.tables[0][3][3], poisson_sum.tables[0][4][1], poisson_sum.tables[0][4][3], poisson_sum.tables[0][5][1], poisson_sum.tables[0][5][3], poisson_sum.tables[0][6][3], poisson_res.params['const'], poisson_res.pvalues['const'],poisson_res.params['step'],poisson_res.pvalues['step'],poisson_res.params['step2'],poisson_res.pvalues['step2'],poisson_res.params['step3'],poisson_res.pvalues['step3'],None,None],index=modelOutputCols)
outputDF = outputDF.append(row, ignore_index=True)
row = pd.Series(["NegBin3",radCutOff, nbin_sum.tables[0][0][1], nbin_sum.tables[0][0][3], nbin_sum.tables[0][1][1], nbin_sum.tables[0][1][3], nbin_sum.tables[0][2][1], nbin_sum.tables[0][2][3], nbin_sum.tables[0][3][1], nbin_sum.tables[0][3][3], nbin_sum.tables[0][4][1], nbin_sum.tables[0][4][3], nbin_sum.tables[0][5][1], nbin_sum.tables[0][5][3], nbin_sum.tables[0][6][3], nbin_res.params['const'],nbin_res.pvalues['const'],nbin_res.params['step'],nbin_res.pvalues['step'],nbin_res.params['step2'],nbin_res.pvalues['step2'],nbin_res.params['step3'],nbin_res.pvalues['step3'],nbin_res.params['alpha'],nbin_res.pvalues['alpha']],index=modelOutputCols)
outputDF = outputDF.append(row, ignore_index=True)
return outputDF
modelOutputCols=['RegName','RadCutOff', 'Dep Var','No. Obs','Model','DF Res','Method','DF Model','Date','PseudoRsq','Time','LL','Converge','LLNull','LLR Pval','CoefCons','PvalCons','CoefStep','PvalStep','CoefStep2','PvalStep2','CoefStep3','PvalStep3','Alpha','PvalAlpha']
modelOutput = pd.DataFrame(columns=modelOutputCols)
modelOutput = doRegs(2,modelOutput)
modelOutput = doRegs(3,modelOutput)
# Export output to CSV
if saveFigs:
fileName = "{0}_Percolation_Innov_RegsOutput.csv".format(fileTime)
modelOutput.to_csv(fileName, na_rep="None")
Simulation Test 1 - Simple and Complex Technology Space
Test 1 compares different values of the monopoly profit decay rate and patent duration on both a simple and complex technology space. The monopoly profit decay rate captures monopoly power, which determines how quickly the competitive fringe erodes the initial monopoly power of innovators. Patent duration is the length of patent protection. In the first test, Test 1.a, the mean of the log normal resistance distribution is low, meaning the technology space is relatively easy to navigate.
The results for Test 1.a suggest that when the technology space is simple and monopoly power is relatively low, the stifling effects of patents blocking off paths of innovative search are outweighed by the benefits of additional resources those monopoly rents provide. As monopoly power increases, on the other hand, firms can effectively navigate the technology space without the additional resources granted by patents.
In [47]:
# Generate list of piDecay values for the test (with repeated values)
piDecay = test1A[(test1A.patentLife==testPatentLife[0])].piDecayRate
# Find maxBPF for each simulation
finalMaxBPFNoPat = test1A[(test1A.doPatent==False)].maxBPF
final5MaxBPF = test1A[(test1A.patentLife==testPatentLife[0]) & (test1A.doPatent==True)].maxBPF
final200MaxBPF = test1A[(test1A.patentLife==testPatentLife[1]) & (test1A.doPatent==True)].maxBPF
title1 = 'Max BPF Height, Simple Tech Space: \n With and Without Patents, Varied Patent Duration'
title2 = 'Max BPF Height, Simple Tech Space: \n Percent Difference Patent and No Patents'
xlab = 'Monopoly Power'
fig = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, piDecay, title1, title2, xlab, True)
if saveFigs:
figName = "{0}_Test1A_MaxBPF.png".format(fileTime)
plt.savefig(figName)
plt.show()
Test 1.b uses the same parameters as Test 1.a but with a more difficult technology space (higher mean resistance). As the technology space becomes more difficult, innovative performance again increases as monopoly power increases. The absolute levels of BPF height attained in Test 1.b are much lower than those seen Test 1.a, primarily because firms must spend more resources to navigate the technology space.
These results show that when the technology space is complex patents improve innovative performance but in a relatively noisy manner.
In [48]:
piDecay = test1B[(test1B.patentLife==testPatentLife[0])].piDecayRate
finalMaxBPFNoPat = test1B[(test1B.doPatent==False)].maxBPF
final5MaxBPF = test1B[(test1B.patentLife==testPatentLife[0]) & (test1B.doPatent==True)].maxBPF
final200MaxBPF = test1B[(test1B.patentLife==testPatentLife[1]) & (test1B.doPatent==True)].maxBPF
title1 = 'Max BPF Height, Complex Tech Space: \n With and Without Patents, Varied Patent Duration'
title2 = 'Max BPF Height Complex Tech Space: \n Percent Difference With and Without Patents'
xlab = 'Monopoly Power'
fig = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, piDecay, title1, title2, xlab, True)
if saveFigs:
figName = "{0}_Test1B_MaxBPF.png".format(fileTime)
plt.savefig(figName)
plt.show()
Simulation Test 2 - Patent Breadth
The second simulation test examines how the model performs with increased patent breadth. The simulation parameters for Tests 2.a and 2.b are identical to those of Test 1.a and Test 1.b but with a higher patent breadth. Increasing patent breadth increases both the potential profits for the innovating firm as well as the stifling effects associated with blocking innovative search.
Comparing to Test 1.a, lower monopoly power no longer implies that patents improve innovative search when patent duration is high. When the technology space is simple, patents are always stifling with a higher patent breadth. The percentage difference plot shows that for all patent durations and levels of monopoly power performance is nearly always higher without patents.
In [49]:
piDecay = test2A[(test2A.patentLife==testPatentLife[0])].piDecayRate
finalMaxBPFNoPat = test2A[(test2A.doPatent==False)].maxBPF
final5MaxBPF = test2A[(test2A.patentLife==testPatentLife[0]) & (test2A.doPatent==True)].maxBPF
final200MaxBPF = test2A[(test2A.patentLife==testPatentLife[1]) & (test2A.doPatent==True)].maxBPF
title1 = 'Max BPF Height, Simple Tech Space: \n With and Without Patents, High Patent Radius'
title2 = 'Max BPF Height Simple Tech Space: \n Percent Difference, High Patent Radius'
xlab = 'Monopoly Power'
fig = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, piDecay, title1, title2, xlab, True)
if saveFigs:
figName = "{0}_Test2A_MaxBPF.png".format(fileTime)
plt.savefig(figName)
plt.show()
The max BPF plots for Test 2.b are shown below. Again, with a more complex technology space, the case without patents loses its dominance. The behavior of the model seems much more erratic with no conclusive pattern shown in the percentage difference plot. This lack of difference may be due to the fact that when the technology space is complex even the smallest differences in the path of innovative search can have large effects on innovative performance.
In [50]:
piDecay = test2B[(test2B.patentLife==testPatentLife[0])].piDecayRate
finalMaxBPFNoPat = test2B[(test2B.doPatent==False)].maxBPF
final5MaxBPF = test2B[(test2B.patentLife==testPatentLife[0]) & (test2B.doPatent==True)].maxBPF
final200MaxBPF = test2B[(test2B.patentLife==testPatentLife[1]) & (test2B.doPatent==True)].maxBPF
title1 = 'Max BPF Height, Complex Tech Space: \n With and Without Patents, High Patent Radius'
title2 = 'Max BPF Height, Complex Tech Space: \n Percent Difference, High Patent Radius'
xlab = 'Monopoly Power'
fig = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, piDecay, title1, title2, xlab, True)
if saveFigs:
figName = "{0}_Test2B_MaxBPF.png".format(fileTime)
plt.savefig(figName)
plt.show()
The figure below shows the difference between the Tests 1 and 2 across patent breadth. The first plot, for the simple technology space, shows that lower patent breadth significantly outperforms the higher patent breadth. XXXXXXXXXXXXXXXXXXXXXXXXX
In [51]:
piDecay = test1A[(test1A.patentLife==testPatentLife[0])].piDecayRate
finalMaxBPFNoPat1S = test1A[(test1A.doPatent==False)].maxBPF
final5MaxBPF1S = test1A[(test1A.patentLife==testPatentLife[0]) & (test1A.doPatent==True)].maxBPF
final200MaxBPF1S = test1A[(test1A.patentLife==testPatentLife[1]) & (test1A.doPatent==True)].maxBPF
finalMaxBPFNoPat3S = test2A[(test2A.doPatent==False)].maxBPF
final5MaxBPF3S = test2A[(test2A.patentLife==testPatentLife[0]) & (test2A.doPatent==True)].maxBPF
final200MaxBPF3S = test2A[(test2A.patentLife==testPatentLife[1]) & (test2A.doPatent==True)].maxBPF
finalMaxBPFNoPat1C = test1B[(test1B.doPatent==False)].maxBPF
final5MaxBPF1C = test1B[(test1B.patentLife==testPatentLife[0]) & (test1B.doPatent==True)].maxBPF
final200MaxBPF1C = test1B[(test1B.patentLife==testPatentLife[1]) & (test1B.doPatent==True)].maxBPF
finalMaxBPFNoPat3C = test2B[(test2B.doPatent==False)].maxBPF
final5MaxBPF3C = test2B[(test2B.patentLife==testPatentLife[0]) & (test2B.doPatent==True)].maxBPF
final200MaxBPF3C = test2B[(test2B.patentLife==testPatentLife[1]) & (test2B.doPatent==True)].maxBPF
# Define fitted lines
pf1 = np.polyfit(piDecay, finalMaxBPFNoPat1S, 4)
p1 = np.poly1d(pf1)
pf2 = np.polyfit(piDecay, final5MaxBPF1S,4)
p2 = np.poly1d(pf2)
pf4 = np.polyfit(piDecay,final200MaxBPF1S,4)
p4 = np.poly1d(pf4)
pf5 = np.polyfit(piDecay, finalMaxBPFNoPat3S, 4)
p5 = np.poly1d(pf5)
pf6 = np.polyfit(piDecay,final5MaxBPF3S,4)
p6 = np.poly1d(pf6)
pf8 = np.polyfit(piDecay,final200MaxBPF3S,4)
p8 = np.poly1d(pf8)
pf9 = np.polyfit(piDecay, finalMaxBPFNoPat1C, 4)
p9 = np.poly1d(pf9)
pf10 = np.polyfit(piDecay,final5MaxBPF1C,4)
p10 = np.poly1d(pf10)
pf12 = np.polyfit(piDecay,final200MaxBPF1C,4)
p12 = np.poly1d(pf12)
pf13 = np.polyfit(piDecay, finalMaxBPFNoPat3C, 4)
p13 = np.poly1d(pf13)
pf14 = np.polyfit(piDecay,final5MaxBPF3C,4)
p14 = np.poly1d(pf14)
pf16 = np.polyfit(piDecay,final200MaxBPF3C,4)
p16 = np.poly1d(pf16)
# Create x axis variable
x2 = np.linspace(round(piDecay.values[0],3), round(piDecay.values[-1],3), 10)
# Create plots
fig = plt.figure(figsize=(13,5))
ax1 = fig.add_subplot(121)
ax1.plot(x2, (p2(x2)-p6(x2))/((p2(x2)+p6(x2))/2), lw=2, ls='--', c='r')
ax1.plot(x2, (p4(x2)-p8(x2))/((p4(x2)+p8(x2))/2), lw=2, ls=':', c='cyan')
ax1.axhline(0, color='black')
ax1.invert_xaxis()
ax1.set_title('Percent Difference Low-High Patent Radius: \n Simple Technology Space')
ax1.set_ylabel('Percent Diff (PatRad1-PatRad3)')
ax1.set_xlabel('Monopoly Power')
labels = ['Low Duration', 'High Duration']
ax1.legend(labels, loc=0)
ax2 = fig.add_subplot(122)
ax2.plot(x2, (p10(x2)-p14(x2))/((p10(x2)+p14(x2))/2), lw=2, ls='--', c='r')
ax2.plot(x2, (p12(x2)-p16(x2))/((p12(x2)+p16(x2))/2), lw=2, ls=':', c='cyan')
ax2.axhline(0, color='black')
ax2.invert_xaxis()
ax2.set_title('Percent Difference Low-High Patent Radius: \n Complex Technology Space')
ax2.set_ylabel('Percent Diff (PatRad1-PatRad3)')
ax2.set_xlabel('Monopoly Power')
labels = ['Low Duration', 'High Duration']
ax2.legend(labels, loc=0)
if saveFigs:
figName = "{0}_Test1_Test2_Diff.png".format(fileTime)
plt.savefig(figName)
plt.show()
Simulation Test 3 - Technology Space Complexity
Tests 3.a and 3.b more closely examine how the complexity of the technology space affects innovative performance. Adjusting the mean of the log-normal distribution for cell resistance makes the technology space relatively more or less difficult to navigate. Test 3.a simulates with strong monopoly power and Test 3.b with weak monopoly power. The figure below shows the maximum BPF reached at the end of each run across different values of the mean resistance for low and high patent durations. The results for Test 3.a suggest that with strong monopoly power, as technological complexity increases, innovative performance with patents increases. There is a critical level of complexity at which a decrease in complexity leads to the no-patent case doing better and an increase in complexity leads to the patent case doing better. For the conditions of Test 3.a, the stifling effects of patents outweigh the benefits of increased monopoly rents up to a certain point of technological complexity.
In [52]:
resistMu = test3A[(test3A.patentLife==testPatentLife[0])].resistMu
# Decay test max BPF plots
finalMaxBPFNoPat = test3A[(test3A.doPatent==False)].maxBPF
final5MaxBPF = test3A[(test3A.patentLife==testPatentLife[0]) & (test3A.doPatent==True)].maxBPF
final200MaxBPF = test3A[(test3A.patentLife==testPatentLife[1]) & (test3A.doPatent==True)].maxBPF
title1 = 'Max BPF Height, High Monopoly Power: \n With and Without Patents, Varied Mean Resistance'
title2 = 'Max BPF Height, High Monopoly Power: \n Percent Difference With and Without Patents'
xlab = 'Mean of Log Normal Resistance Dist.'
fig = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, resistMu, title1, title2, xlab, False)
if saveFigs:
figName = "{0}_Test3A_MaxBPF.png".format(fileTime)
plt.savefig(figName)
plt.show()
Test 3.b shows the opposite result when monopoly power is weak. Shown in the figure below, when monopoly power is relatively low the simulation with patents produces more innovation. Again, the gap in performance between the two cases decreases as the technology space becomes more complex. In contrast to Test 3.a, when monopoly power is low, increased patent duration improves innovative search across a wide range of technological complexities. With weaker monopoly power, the additional resources are necessary to propel firms through the technology space, with the benefits of those additional resources outweighing any of the negative effects of blocking.
In [53]:
resistMu = test3B[(test3B.patentLife==testPatentLife[0])].resistMu
# Decay test max BPF plots
finalMaxBPFNoPat = test3B[(test3B.doPatent==False)].maxBPF
final5MaxBPF = test3B[(test3B.patentLife==testPatentLife[0]) & (test3B.doPatent==True)].maxBPF
final200MaxBPF = test3B[(test3B.patentLife==testPatentLife[1]) & (test3B.doPatent==True)].maxBPF
title1 = 'Max BPF Height, Low Monopoly Power: \n With and Without Patents, Varied Mean Resistance'
title2 = 'Max BPF Height, Low Monopoly Power: \n Percent Difference With and Without Patents'
xlab = 'Mean of Log Normal Resistance Dist.'
fig = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, resistMu, title1, title2, xlab, False)
if saveFigs:
figName = "{0}_Test3B_MaxBPF.png".format(fileTime)
plt.savefig(figName)
plt.show()
Simulation Test 4 - Technology Space Complexity
Tests 4.a and 4.b show how changes in the size of monopoly rents affect innovative activity. Higher monopoly rents imply that firms will have more resources to invest in R&D, both with and without patents. Test 4.a uses a simple technology space, and Test 4.b uses a complex technology space. Similar to the results from Test 1.a and Test 2.a, the figure shows that when the technology space is simple, the case without patents generates more innovative activity relative to the case with patents as monopoly power increases. The percent difference between the case without and with patents for each level of monopoly rents is shown below. The differences are positive for higher values of monopoly power, and the gap is increasing in monopoly rents.
In [54]:
testMonopProfit = [1, 3, 5]
piDecay = test4A[(test4A.monopProfit==testMonopProfit[0]) & (test4A.doPatent==False)].piDecayRate
finalMaxBPFNoPat1 = test4A[(test4A.doPatent==False) & (test4A.monopProfit==testMonopProfit[0])].maxBPF
finalMaxBPFNoPat3 = test4A[(test4A.doPatent==False) & (test4A.monopProfit==testMonopProfit[1])].maxBPF
finalMaxBPFNoPat5 = test4A[(test4A.doPatent==False) & (test4A.monopProfit==testMonopProfit[2])].maxBPF
final1MaxBPF = test4A[(test4A.monopProfit==testMonopProfit[0]) & (test4A.doPatent==True)].maxBPF
final3MaxBPF = test4A[(test4A.monopProfit==testMonopProfit[1]) & (test4A.doPatent==True)].maxBPF
final5MaxBPF = test4A[(test4A.monopProfit==testMonopProfit[2]) & (test4A.doPatent==True)].maxBPF
pf1 = np.polyfit(piDecay, finalMaxBPFNoPat1, 4)
p1 = np.poly1d(pf1)
pf2 = np.polyfit(piDecay, finalMaxBPFNoPat3, 4)
p2 = np.poly1d(pf2)
pf3 = np.polyfit(piDecay, finalMaxBPFNoPat5, 4)
p3 = np.poly1d(pf3)
pf4 = np.polyfit(piDecay,final1MaxBPF,4)
p4 = np.poly1d(pf4)
pf5 = np.polyfit(piDecay, final3MaxBPF, 4)
p5 = np.poly1d(pf5)
pf6 = np.polyfit(piDecay,final5MaxBPF,4)
p6 = np.poly1d(pf6)
x2 = np.linspace(0.01, 0.1, 10)
fig = plt.figure(figsize=(13,5))
ax1 = fig.add_subplot(121)
ax1.plot(x2, p1(x2), lw=2, c='k')
ax1.plot(x2, p2(x2), lw=2, ls='--', c='m')
ax1.plot(x2, p3(x2), lw=2, ls='-.', c='orange')
ax1.invert_xaxis()
ax1.set_title('Max BPF Height Without Patents: \n Simple Tech Space, Varied Monopoly Rents')
ax1.set_ylabel('Max BPF Height')
ax1.set_xlabel('Monopoly Power')
labels = ['Low Rents', 'Moderate Rents', 'High Rents']
ax1.legend(labels, loc=0)
ax2 = fig.add_subplot(122)
ax2.plot(x2, p4(x2), lw=2, c='k')
ax2.plot(x2, p5(x2), lw=2, ls='--', c='m')
ax2.plot(x2, p6(x2), lw=2, ls='-.', c='orange')
ax2.invert_xaxis()
ax2.set_title('Max BPF Height With Patents: \n Simple Tech Space, Varied Monopoly Rents')
ax2.set_ylabel('Max BPF Height')
ax2.set_xlabel('Monopoly Power')
ax2.set_ylim((0,6000))
labels = ['Low Rents', 'Moderate Rents', 'High Rents']
ax2.legend(labels, loc=0)
if saveFigs:
figName = "{0}_Test4A_MaxBPF.png".format(fileTime)
plt.savefig(figName)
plt.show()
fig = plt.figure(figsize=(6,5))
ax1 = fig.add_subplot(111)
ax1.plot(x2, (p1(x2)-p4(x2))/((p1(x2)+p4(x2))/2), lw=2, c='k')
ax1.plot(x2, (p2(x2)-p5(x2))/((p2(x2)+p5(x2))/2), lw=2, ls='--', c='m')
ax1.plot(x2, (p3(x2)-p6(x2))/((p3(x2)+p6(x2))/2), lw=2, ls='-.', c='orange')
ax1.invert_xaxis()
ax1.axhline(0, color='black')
ax1.set_title('Max BPF Height, Simple Tech Space: \n Percent Difference With and Without Patents')
ax1.set_ylabel('Max BPF Height')
ax1.set_xlabel('Monopoly Power')
labels = ['Low Rents', 'Moderate Rents', 'High Rents']
ax1.legend(labels, loc=0)
if saveFigs:
figName = "{0}_Test4A_Diff.png".format(fileTime)
plt.savefig(figName)
plt.show()
Test 4.b runs simulations identical to those in Test 4.a but with a complex technology space. In the plot below we see a much more drastic improvement in model performance for the highest values of monopoly power. The figure below shows the percent difference between Test 4.b runs without patents and with patents. Unlike in Test 4.a, with a more complex technology space the percent differences are much more volatile. The vast dispersion in outcomes can be in part attributed to the increased importance of path dependence, a result familiar to historians of innovation and technological change. When the technology space is complex a small change in direction can drastically improve innovative outcomes.
In [55]:
piDecay = test4B[(test4B.monopProfit==testMonopProfit[0]) & (test4B.doPatent==False)].piDecayRate
finalMaxBPFNoPat1 = test4B[(test4B.doPatent==False) & (test4B.monopProfit==testMonopProfit[0])].maxBPF
finalMaxBPFNoPat3 = test4B[(test4B.doPatent==False) & (test4B.monopProfit==testMonopProfit[1])].maxBPF
finalMaxBPFNoPat5 = test4B[(test4B.doPatent==False) & (test4B.monopProfit==testMonopProfit[2])].maxBPF
final1MaxBPF = test4B[(test4B.monopProfit==testMonopProfit[0]) & (test4B.doPatent==True)].maxBPF
final3MaxBPF = test4B[(test4B.monopProfit==testMonopProfit[1]) & (test4B.doPatent==True)].maxBPF
final5MaxBPF = test4B[(test4B.monopProfit==testMonopProfit[2]) & (test4B.doPatent==True)].maxBPF
pf1 = np.polyfit(piDecay, finalMaxBPFNoPat1, 3)
p1 = np.poly1d(pf1)
pf2 = np.polyfit(piDecay, finalMaxBPFNoPat3, 3)
p2 = np.poly1d(pf2)
pf3 = np.polyfit(piDecay, finalMaxBPFNoPat5, 3)
p3 = np.poly1d(pf3)
pf4 = np.polyfit(piDecay,final1MaxBPF, 3)
p4 = np.poly1d(pf4)
pf5 = np.polyfit(piDecay, final3MaxBPF, 3)
p5 = np.poly1d(pf5)
pf6 = np.polyfit(piDecay,final5MaxBPF, 3)
p6 = np.poly1d(pf6)
x2 = np.linspace(0.01, 0.1, 10)
fig = plt.figure(figsize=(13,5))
ax1 = fig.add_subplot(121)
ax1.plot(x2, p1(x2), lw=2, c='k')
ax1.plot(x2, p2(x2), lw=2, ls='--', c='m')
ax1.plot(x2, p3(x2), lw=2, ls='-.', c='orange')
ax1.invert_xaxis()
ax1.set_title('Max BPF Height Without Patents: \n Complex Tech Space, Varied Monopoly Rents')
ax1.set_ylabel('Max BPF Height')
ax1.set_xlabel('Monopoly Power')
labels = ['Low Rents', 'Moderate Rents', 'High Rents']
ax1.legend(labels, loc=0)
ax2 = fig.add_subplot(122)
ax2.plot(x2, p4(x2), lw=2, c='k')
ax2.plot(x2, p5(x2), lw=2, ls='--', c='m')
ax2.plot(x2, p6(x2), lw=2, ls='-.', c='orange')
ax2.invert_xaxis()
ax2.set_title('Max BPF Height With Patents: \n Complex Tech Space, Varied Monopoly Rents')
ax2.set_ylabel('Max BPF Height')
ax2.set_xlabel('Monopoly Power')
ax2.set_ylim((0,2000))
labels = ['Low Rents', 'Moderate Rents', 'High Rents']
ax2.legend(labels, loc=0)
if saveFigs:
figName = "{0}_Test4B_MaxBPF.png".format(fileTime)
plt.savefig(figName)
plt.show()
fig = plt.figure(figsize=(6,5))
ax1 = fig.add_subplot(111)
ax1.plot(x2, (p1(x2)-p4(x2))/((p1(x2)+p4(x2))/2), lw=2, c='k')
ax1.plot(x2, (p2(x2)-p5(x2))/((p2(x2)+p5(x2))/2), lw=2, ls='--', c='m')
ax1.plot(x2, (p3(x2)-p6(x2))/((p3(x2)+p6(x2))/2), lw=2, ls='-.', c='orange')
ax1.invert_xaxis()
ax1.axhline(0, color='black')
ax1.set_title('Max BPF Height, Complex Tech Space: \n Percent Difference With and Without Patents')
ax1.set_ylabel('Max BPF Height')
ax1.set_xlabel('Monopoly Power')
labels = ['Low Rents', 'Moderate Rents', 'High Rents']
ax1.legend(labels, loc=0)
if saveFigs:
figName = "{0}_Test4B_Diff.png".format(fileTime)
plt.savefig(figName)
plt.show()