This ipython notebook produces the graphs used to analyze the output of an agent-based model of the innovation process. File names below need to be modified to point to the output files to be analyzed, order in the list matters as they are picked up by the pandas data frames.
In [67]:
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 [68]:
# Lists of file names for simulation results csvs
testdatafiles = ['2014-10-20 0551_39_test1a_percol_rundata.csv',
'2014-10-20 0551_39_test1b_percol_rundata.csv',
'2014-10-19 2058_17_test2a_percol_rundata.csv',
'2014-10-19 2058_17_test2b_percol_rundata.csv',
'2014-10-23 1330_54_test3a_percol_rundata.csv',
'2014-10-23 1330_54_test3b_percol_rundata.csv',
'2014-10-23 1331_26_test4a_percol_rundata.csv',
'2014-10-23 1331_26_test4b_percol_rundata.csv',
'2014-10-26 2039_34_microSim_percol_rundata.csv',
'2015-01-11 2141_56_sim_typical_0_percol_innovdata.csv',
'2015-01-11 2141_56_sim_typical_0_percol_stepdata.csv',
]
testparamfiles = ['2014-10-20 0551_39_test1a_percol_params.csv',
'2014-10-20 0551_39_test1b_percol_params.csv',
'2014-10-19 2058_17_test2a_percol_params.csv',
'2014-10-19 2058_17_test2b_percol_params.csv',
'2014-10-23 1330_54_test3a_percol_params.csv',
'2014-10-23 1330_54_test3b_percol_params.csv',
'2014-10-23 1331_26_test4a_percol_params.csv',
'2014-10-23 1331_26_test4b_percol_params.csv',
'2014-10-26 2039_34_microSim_percol_params.csv',
'2015-01-11 2141_56_sim_typical_0_percol_params.csv'
]
test1A = pd.merge(pd.DataFrame.from_csv(testdatafiles[0],index_col=False),pd.DataFrame.from_csv(testparamfiles[0],index_col=False), on='run')
test1B = pd.merge(pd.DataFrame.from_csv(testdatafiles[1],index_col=False),pd.DataFrame.from_csv(testparamfiles[1],index_col=False), on='run')
test2A = pd.merge(pd.DataFrame.from_csv(testdatafiles[2],index_col=False),pd.DataFrame.from_csv(testparamfiles[2],index_col=False), on='run')
test2B = pd.merge(pd.DataFrame.from_csv(testdatafiles[3],index_col=False),pd.DataFrame.from_csv(testparamfiles[3],index_col=False), on='run')
test3A = pd.merge(pd.DataFrame.from_csv(testdatafiles[4],index_col=False),pd.DataFrame.from_csv(testparamfiles[4],index_col=False), on='run')
test3B = pd.merge(pd.DataFrame.from_csv(testdatafiles[5],index_col=False),pd.DataFrame.from_csv(testparamfiles[5],index_col=False), on='run')
test4A = pd.merge(pd.DataFrame.from_csv(testdatafiles[6],index_col=False),pd.DataFrame.from_csv(testparamfiles[6],index_col=False), on='run')
test4B = pd.merge(pd.DataFrame.from_csv(testdatafiles[7],index_col=False),pd.DataFrame.from_csv(testparamfiles[7],index_col=False), on='run')
microSim = pd.merge(pd.DataFrame.from_csv(testdatafiles[8],index_col=False),pd.DataFrame.from_csv(testparamfiles[8],index_col=False), on='run')
invDF = pd.DataFrame.from_csv(testdatafiles[9],index_col=False)
stepDF = pd.DataFrame.from_csv(testdatafiles[10],index_col=False)
microSim['lmaxBPF']=np.log(microSim['maxBPF'])
fileTime = time.strftime("%Y-%m-%d %H%M_%S")
testPatentLife = [5, 200]
testMonopProfit = [1, 3, 5]
radCutOff = 2
numSteps = 3000
In [69]:
# 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
fig1 = plt.figure(figsize=(5,4))
ax1 = fig1.add_subplot(111)
ax1.plot(x2, pn(x2), lw=2.5, c='black')
ax1.plot(x2, p5(x2), lw=2.5, ls='--', c='blue')
ax1.plot(x2, p200(x2), lw=2.5, ls=':', c='red')
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)
fig2 = plt.figure(figsize=(5,4))
ax2 = fig2.add_subplot(111)
ax2.plot(x2, (p5(x2)-pn(x2))/((pn(x2)+p5(x2))/2), lw=2, ls='--', c='blue')
ax2.plot(x2, (p200(x2)-pn(x2))/((pn(x2)+p200(x2))/2), lw=2, ls=':', c='red')
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 fig1, fig2
In [70]:
# 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 [71]:
# Select the innovation sizes for one of the runs
innovSizes = invDF.innovSize
# Plot histogram of innovation size
fig = plt.figure(figsize=(5,4))
ax1 = fig.add_subplot(111)
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)
if saveFigs:
figName = "TypicalRun_InnovHist.png"
plt.savefig(figName)
plt.show()
# Plot log series
fig = plt.figure(figsize=(5,4))
ax2 = fig.add_subplot(111)
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 = "TypicalRun_logplot.png"
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 [72]:
# Select BPF Change time series for one of the runs
BPFChange = stepDF[stepDF.run==0].BPFChange
# Plot BPF Change time series
fig = plt.figure(figsize=(5,4))
ax1 = fig.add_subplot(111)
ax1.plot(BPFChange)
plt.xlabel('Step')
plt.ylabel('BPF Change')
plt.title('BPF Change Time Series')
if saveFigs:
figName = "TypicalRun_BPFChange.png".format(fileTime)
plt.savefig(figName)
plt.show()
# BPF Change lag plot
fig = plt.figure(figsize=(5,4))
ax1 = fig.add_subplot(111)
lag_plot(BPFChange)
plt.xlabel('BPFChange(t)')
plt.ylabel('BPFChange(t + 1)')
plt.title('BPF Change Lag Plot')
if saveFigs:
figName = "TypicalRun_BPFlagplot.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 [73]:
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 = invDF[(invDF.run==0) & (invDF.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.
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 [74]:
# 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'
fig1, fig2 = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, piDecay, title1, title2, xlab, True)
if saveFigs:
fig1.savefig("Test1A_MaxBPF.png",bbox_inches='tight')
fig2.savefig("Test1A_DiffBPF.png",bbox_inches='tight')
In [75]:
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'
fig1, fig2 = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, piDecay, title1, title2, xlab, True)
if saveFigs:
fig1.savefig("Test1B_MaxBPF.png",bbox_inches='tight')
fig2.savefig("Test1B_DiffBPF.png",bbox_inches='tight')
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 [76]:
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'
fig1, fig2 = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, piDecay, title1, title2, xlab, True)
if saveFigs:
fig1.savefig("Test2A_MaxBPF.png",bbox_inches='tight')
fig2.savefig("Test2A_DiffBPF.png",bbox_inches='tight')
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 [77]:
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'
fig1, fig2 = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, piDecay, title1, title2, xlab, True)
if saveFigs:
fig1.savefig("Test2B_MaxBPF.png",bbox_inches='tight')
fig2.savefig("Test2B_DiffBPF.png",bbox_inches='tight')
In [78]:
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=(5,4))
ax1 = fig.add_subplot(111)
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)
if saveFigs:
plt.savefig("Test1_Test2_Sim_Diff.png",bbox_inches='tight')
plt.show()
fig = plt.figure(figsize=(5,4))
ax2 = fig.add_subplot(111)
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:
plt.savefig("Test1_Test2_Cmplx_Diff.png",bbox_inches='tight')
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 [79]:
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.'
fig1,fig2 = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, resistMu, title1, title2, xlab, False)
if saveFigs:
fig1.savefig("Test3A_MaxBPF.png",bbox_inches='tight')
fig2.savefig("Test3A_DiffBPF.png",bbox_inches='tight')
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 [80]:
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.'
fig1, fig2 = doMaxBPFPlot(finalMaxBPFNoPat, final5MaxBPF, final200MaxBPF, resistMu, title1, title2, xlab, False)
if saveFigs:
fig1.savefig("Test3B_MaxBPF.png",bbox_inches='tight')
fig2.savefig("Test3B_DiffBPF.png",bbox_inches='tight')
Micro-Simulations
In addition to the structured simulation tests described in the previous sections, a large number of micro-simulations were run using parameters drawn from uniform distributions. The range of the distribution for each parameter are similar to those described in the three simulation tests in the previous sections. The Figure below shows scatter plots of the maximum BPF reached against the complexity of the technology space and the left panel shows the maximum BPF reached against the level of firms' monopoly power.
In [81]:
fig = plt.figure(figsize=(5,4))
ax1 = fig.add_subplot(111)
ax1.scatter(microSim[microSim.doPatent==True]['resistMu'],microSim[microSim.doPatent==True]['maxBPF'],color='b',s=10,alpha=0.5)
ax1.scatter(microSim[microSim.doPatent==False]['resistMu'],microSim[microSim.doPatent==False]['maxBPF'],color='r',s=10,alpha=0.5,marker='^')
ax1.set_ylabel('Max BPF Height')
ax1.set_xlabel('Mean of Log Normal Resistance Dist.')
ax1.set_xlim(2.0,4.0)
ax1.set_ylim(0,500)
ax1.set_title('Micro-Simulation Max BPF Height vs. Mean Resistance')
labels = ['With Patents', 'Without Patents']
ax1.legend(labels)
if saveFigs:
figName = "microsim_maxbpf_resistmu.png".format(fileTime)
plt.savefig(figName)
plt.show()
fig = plt.figure(figsize=(5,4))
ax1 = fig.add_subplot(111)
ax1.scatter(microSim[microSim.doPatent==True]['piDecayRate'],microSim[microSim.doPatent==True]['maxBPF'],color='b',s=10,alpha=0.5)
ax1.scatter(microSim[microSim.doPatent==False]['piDecayRate'],microSim[microSim.doPatent==False]['maxBPF'],color='r',s=10,alpha=0.5,marker='^')
ax1.set_ylabel('Max BPF Height')
ax1.set_xlabel('Monopoly Power')
ax1.set_xlim(0.01,0.1)
ax1.set_ylim(0,500)
ax1.invert_xaxis()
ax1.set_title('Micro-Simulation Max BPF Height vs. Monopoly Power')
labels = ['With Patents', 'Without Patents']
ax1.legend(labels, loc='upper left')
if saveFigs:
figName = "microsim_maxbpf_pidecay.png"
plt.savefig(figName)
plt.show()
.
Closer examination of the distributions of innovative performance in simulations with and without patents reveals the conditions under which patents can deter innovative activity. We will focus on the two parameters of interest, the complexity of the technology space and monopoly power. The figures below shows the distribution of logged maximum BPF height reached in simulations with and without patents and a simple technology space.
In [82]:
labels = ['Patents','No Patents']
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(3,4))
dataPats = microSim[(microSim.doPatent==True)]['maxBPF']
dataNoPats = microSim[(microSim.doPatent==False)]['maxBPF']
data = np.array([dataPats,dataNoPats])
axes.boxplot(data,labels=labels,showmeans=True)
axes.set_ylabel('Max BPF')
axes.set_title('Max BPF Height\nOverall')
plt.show()
dataPats = microSim[(microSim.doPatent==True)&(microSim.resistMu<=2.5)]['maxBPF']
dataNoPats = microSim[(microSim.doPatent==False)&(microSim.resistMu<=2.5)]['maxBPF']
data = np.array([dataPats,dataNoPats])
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(3,4))
axes.boxplot(data,labels=labels,showmeans=True)
axes.set_title('Max BPF Height\nSimple Technology Space')
axes.set_ylabel('Max BPF')
if saveFigs:
figName = "microsim_bplot_simtech.png"
plt.savefig(figName,bbox_inches='tight')
plt.show()
dataPats = microSim[(microSim.doPatent==True)&(microSim.resistMu>=3.5)]['maxBPF']
dataNoPats = microSim[(microSim.doPatent==False)&(microSim.resistMu>=3.5)]['maxBPF']
data = np.array([dataPats,dataNoPats])
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(3,4))
axes.boxplot(data,labels=labels,showmeans=True)
axes.set_title('Max BPF Height\nComplex Technology Space')
axes.set_ylabel('Max BPF')
if saveFigs:
figName = "microsim_bplot_cmplxtech.png"
plt.savefig(figName,bbox_inches='tight')
plt.show()
dataPats = microSim[(microSim.doPatent==True)&(microSim.piDecayRate<=0.03)]['maxBPF']
dataNoPats = microSim[(microSim.doPatent==False)&(microSim.piDecayRate<=0.03)]['maxBPF']
data = np.array([dataPats,dataNoPats])
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(3,4))
axes.boxplot(data,labels=labels,showmeans=True)
axes.set_title('Max BPF Height\nHigh Monopoly Power')
axes.set_ylabel('Max BPF')
if saveFigs:
figName = "microsim_bplot_himonop.png"
plt.savefig(figName,bbox_inches='tight')
plt.show()
dataPats = microSim[(microSim.doPatent==True)&(microSim.piDecayRate>=0.08)]['maxBPF']
dataNoPats = microSim[(microSim.doPatent==False)&(microSim.piDecayRate>=0.08)]['maxBPF']
data = np.array([dataPats,dataNoPats])
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(3,4))
axes.boxplot(data,labels=labels,showmeans=True)
axes.set_title('Max BPF Height\nLow Monopoly Power')
axes.set_ylabel('Max BPF')
if saveFigs:
figName = "microsim_bplot_lowmonop.png"
plt.savefig(figName,bbox_inches='tight')
plt.show()
In [83]:
labels=['Patents','No Patents']
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(5,4))
microSim[(microSim.doPatent==True)&(microSim.resistMu<=2.5)]['lmaxBPF'].plot(kind='kde',ax=axes,c='blue',lw=2)
patMean = mean(microSim[(microSim.doPatent==True)&(microSim.resistMu<=2.5)]['lmaxBPF'])
microSim[(microSim.doPatent==False)&(microSim.resistMu<=2.5)]['lmaxBPF'].plot(kind='kde',ax=axes,c='green',lw=2)
noPatMean = mean(microSim[(microSim.doPatent==False)&(microSim.resistMu<=2.5)]['lmaxBPF'])
axes.set_title('Distribution of Log Max BPF\nSimple Technology Space')
axes.set_ylabel('Density')
axes.set_xlabel('Log Max BPF')
axes.legend(labels)
axes.axvline(x=patMean,c='blue',ls='--')
axes.axvline(x=noPatMean,c='green',ls='--')
axes.grid(False)
if saveFigs:
figName = "microsim_dist_simtech.png"
plt.savefig(figName,bbox_inches='tight')
plt.show()
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(5,4))
microSim[(microSim.doPatent==True)&(microSim.resistMu>=3.5)]['lmaxBPF'].plot(kind='kde',ax=axes,c='blue',lw=2)
patMean = mean(microSim[(microSim.doPatent==True)&(microSim.resistMu>=3.5)]['lmaxBPF'])
microSim[(microSim.doPatent==False)&(microSim.resistMu>=3.5)]['lmaxBPF'].plot(kind='kde',ax=axes,c='green',lw=2)
noPatMean = mean(microSim[(microSim.doPatent==False)&(microSim.resistMu>=3.5)]['lmaxBPF'])
axes.set_title('Distribution of Log Max BPF\nComplex Technology Space')
axes.set_ylabel('Density')
axes.set_xlabel('Log Max BPF')
axes.legend(labels)
axes.axvline(x=patMean,c='blue',ls='--')
axes.axvline(x=noPatMean,c='green',ls='--')
axes.grid(False)
if saveFigs:
figName = "microsim_dist_cmplxtech.png"
plt.savefig(figName,bbox_inches='tight')
plt.show()
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(5,4))
microSim[(microSim.doPatent==True)&(microSim.piDecayRate<=0.03)]['lmaxBPF'].plot(kind='kde',ax=axes,c='blue',lw=2)
patMean= mean(microSim[(microSim.doPatent==True)&(microSim.piDecayRate<=0.03)]['lmaxBPF'])
microSim[(microSim.doPatent==False)&(microSim.piDecayRate<=0.03)]['lmaxBPF'].plot(kind='kde',ax=axes,c='green',lw=2)
noPatMean= mean(microSim[(microSim.doPatent==False)&(microSim.piDecayRate<=0.03)]['lmaxBPF'])
axes.set_title('Distribution of Log Max BPF\nHigh Monopoly Power')
axes.set_ylabel('Density')
axes.set_xlabel('Log Max BPF')
axes.legend(labels)
axes.axvline(x=patMean,c='blue',ls='--')
axes.axvline(x=noPatMean,c='green',ls='--')
axes.grid(False)
if saveFigs:
figName = "microsim_dist_himonop.png"
plt.savefig(figName,bbox_inches='tight')
plt.show()
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(5,4))
microSim[(microSim.doPatent==True)&(microSim.piDecayRate>=0.08)]['lmaxBPF'].plot(kind='kde',ax=axes,c='blue',lw=2)
patMean = mean(microSim[(microSim.doPatent==True)&(microSim.piDecayRate>=0.08)]['lmaxBPF'])
microSim[(microSim.doPatent==False)&(microSim.piDecayRate>=0.08)]['lmaxBPF'].plot(kind='kde',ax=axes,c='green',lw=2)
noPatMean = mean(microSim[(microSim.doPatent==False)&(microSim.piDecayRate>=0.08)]['lmaxBPF'])
axes.set_title('Distribution of Log Max BPF\nLow Monopoly Power')
axes.set_ylabel('Density')
axes.set_xlabel('Log Max BPF')
axes.legend(labels,loc='upper left')
axes.axvline(x=patMean,c='blue',ls='--')
axes.axvline(x=noPatMean,c='green',ls='--')
axes.grid(False)
if saveFigs:
figName = "microsim_dist_lowmonop.png"
plt.savefig(figName,bbox_inches='tight')
plt.show()
In [84]:
# calculate mean, std, var, z value for normal dist
print "simple count, no patent", len(microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==False)].maxBPF)
print "simple mean, no patent", microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==False)].maxBPF.mean()
print "simple std, no patent", microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==False)].maxBPF.std()
print "simple count, patent", len(microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==True)].maxBPF)
print "simple mean, patent", microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==True)].maxBPF.mean()
print "simple std, patent", microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==True)].maxBPF.std()
noPatVar = microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==False)].maxBPF.var()
noPatMean = microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==False)].maxBPF.mean()
patVar = microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==True)].maxBPF.var()
patMean = microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==True)].maxBPF.mean()
print "noPat> ==> P(z>", float(patMean-noPatMean)/float(sqrt(noPatVar+patVar)),")","\n"
print "complex count, no patent", len(microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==False)].maxBPF)
print "complex mean, no patent", microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==False)].maxBPF.mean()
print "complex std, no patent", microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==False)].maxBPF.std()
print "complex count, patent", len(microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==True)].maxBPF)
print "complex mean, patent", microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==True)].maxBPF.mean()
print "complex std, patent", microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==True)].maxBPF.std()
noPatVar = microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==False)].maxBPF.var()
noPatMean = microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==False)].maxBPF.mean()
patVar = microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==True)].maxBPF.var()
patMean = microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==True)].maxBPF.mean()
print "noPat> ==> P(z>", float(patMean-noPatMean)/float(sqrt(noPatVar+patVar)),")","\n"
print "high monopoly count, no patent",len(microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==False)].maxBPF)
print "high monopoly mean, no patent", microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==False)].maxBPF.mean()
print "high monopoly std, no patent", microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==False)].maxBPF.std()
print "high monopoly count, patent",len(microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==True)].maxBPF)
print "high monopoly mean, patent", microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==True)].maxBPF.mean()
print "high monopoly std, patent", microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==True)].maxBPF.std()
noPatVar = microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==False)].maxBPF.var()
noPatMean = microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==False)].maxBPF.mean()
patVar = microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==True)].maxBPF.var()
patMean = microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==True)].maxBPF.mean()
print "noPat> ==> P(z>", float(patMean-noPatMean)/float(sqrt(noPatVar+patVar)),")", "\n"
print "low monopoly count, no patent", len(microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==False)])
print "low monopoly mean, no patent", microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==False)].maxBPF.mean()
print "low monopoly std, no patent", microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==False)].maxBPF.std()
print "low monopoly count, patent", len(microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==True)])
print "low monopoly mean, patent", microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==True)].maxBPF.mean()
print "low monopoly std, patent", microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==True)].maxBPF.std()
noPatVar = microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==False)].maxBPF.var()
noPatMean = microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==False)].maxBPF.mean()
patVar = microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==True)].maxBPF.var()
patMean = microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==True)].maxBPF.mean()
print "noPat> ==> P(z>", float(patMean-noPatMean)/float(sqrt(noPatVar+patVar)),")", "\n"
In [85]:
# simulate random draws
noPatMaxBPF = list(log(microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==False)].maxBPF))
patMaxBPF = list(log(microSim[(microSim.resistMu<=2.5)&(microSim.doPatent==True)].maxBPF))
numDraws = 100000
draws = []
for i in range(numDraws):
noPat = np.random.choice(noPatMaxBPF)
pat = np.random.choice(patMaxBPF)
if noPat>pat:
draws.append(1)
else:
draws.append(0)
print 'simple tech, nopat>pat', float(sum(draws))/float(len(draws))
noPatMaxBPF = list(log(microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==False)].maxBPF))
patMaxBPF = list(log(microSim[(microSim.resistMu>=3.5)&(microSim.doPatent==True)].maxBPF))
draws = []
for i in range(numDraws):
noPat = np.random.choice(noPatMaxBPF)
pat = np.random.choice(patMaxBPF)
if noPat>pat:
draws.append(1)
else:
draws.append(0)
print 'complex tech, nopat>pat',float(sum(draws))/float(len(draws))
noPatMaxBPF = list(log(microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==False)].maxBPF))
patMaxBPF = list(log(microSim[(microSim.piDecayRate<=0.03)&(microSim.doPatent==True)].maxBPF))
draws = []
for i in range(numDraws):
noPat = np.random.choice(noPatMaxBPF)
pat = np.random.choice(patMaxBPF)
if noPat>pat:
draws.append(1)
else:
draws.append(0)
print 'high monop, nopat>pat',float(sum(draws))/float(len(draws))
noPatMaxBPF = list(log(microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==False)].maxBPF))
patMaxBPF = list(log(microSim[(microSim.piDecayRate>=0.08)&(microSim.doPatent==True)].maxBPF))
draws = []
for i in range(numDraws):
noPat = np.random.choice(noPatMaxBPF)
pat = np.random.choice(patMaxBPF)
if noPat>pat:
draws.append(1)
else:
draws.append(0)
print 'low monop, nopat>pat',float(sum(draws))/float(len(draws))