Percolating Patents: Balancing the Effects of Patents on Innovation

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 [2]:

    
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

Import Percolation Simulation Data



In [3]:

    
# 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









    



---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-3-7abeb03fae1d> in <module>()
     24                  ]
     25 
---> 26 test1A = pd.merge(pd.DataFrame.from_csv(testdatafiles[0],index_col=False),pd.DataFrame.from_csv(testparamfiles[0],index_col=False), on='run')
     27 test1B = pd.merge(pd.DataFrame.from_csv(testdatafiles[1],index_col=False),pd.DataFrame.from_csv(testparamfiles[1],index_col=False), on='run')
     28 test2A = pd.merge(pd.DataFrame.from_csv(testdatafiles[2],index_col=False),pd.DataFrame.from_csv(testparamfiles[2],index_col=False), on='run')

C:\Users\ngold\AppData\Local\Enthought\Canopy\User\lib\site-packages\pandas\core\frame.pyc in from_csv(cls, path, header, sep, index_col, parse_dates, encoding, tupleize_cols, infer_datetime_format)
   1025                           parse_dates=parse_dates, index_col=index_col,
   1026                           encoding=encoding, tupleize_cols=tupleize_cols,
-> 1027                           infer_datetime_format=infer_datetime_format)
   1028 
   1029     def to_sparse(self, fill_value=None, kind='block'):

C:\Users\ngold\AppData\Local\Enthought\Canopy\User\lib\site-packages\pandas\io\parsers.pyc in parser_f(filepath_or_buffer, sep, dialect, compression, doublequote, escapechar, quotechar, quoting, skipinitialspace, lineterminator, header, index_col, names, prefix, skiprows, skipfooter, skip_footer, na_values, na_fvalues, true_values, false_values, delimiter, converters, dtype, usecols, engine, delim_whitespace, as_recarray, na_filter, compact_ints, use_unsigned, low_memory, buffer_lines, warn_bad_lines, error_bad_lines, keep_default_na, thousands, comment, decimal, parse_dates, keep_date_col, dayfirst, date_parser, memory_map, float_precision, nrows, iterator, chunksize, verbose, encoding, squeeze, mangle_dupe_cols, tupleize_cols, infer_datetime_format, skip_blank_lines)
    463                     skip_blank_lines=skip_blank_lines)
    464 
--> 465         return _read(filepath_or_buffer, kwds)
    466 
    467     parser_f.__name__ = name

C:\Users\ngold\AppData\Local\Enthought\Canopy\User\lib\site-packages\pandas\io\parsers.pyc in _read(filepath_or_buffer, kwds)
    239 
    240     # Create the parser.
--> 241     parser = TextFileReader(filepath_or_buffer, **kwds)
    242 
    243     if (nrows is not None) and (chunksize is not None):

C:\Users\ngold\AppData\Local\Enthought\Canopy\User\lib\site-packages\pandas\io\parsers.pyc in __init__(self, f, engine, **kwds)
    555             self.options['has_index_names'] = kwds['has_index_names']
    556 
--> 557         self._make_engine(self.engine)
    558 
    559     def _get_options_with_defaults(self, engine):

C:\Users\ngold\AppData\Local\Enthought\Canopy\User\lib\site-packages\pandas\io\parsers.pyc in _make_engine(self, engine)
    692     def _make_engine(self, engine='c'):
    693         if engine == 'c':
--> 694             self._engine = CParserWrapper(self.f, **self.options)
    695         else:
    696             if engine == 'python':

C:\Users\ngold\AppData\Local\Enthought\Canopy\User\lib\site-packages\pandas\io\parsers.pyc in __init__(self, src, **kwds)
   1059         kwds['allow_leading_cols'] = self.index_col is not False
   1060 
-> 1061         self._reader = _parser.TextReader(src, **kwds)
   1062 
   1063         # XXX

C:\Users\ngold\AppData\Local\Enthought\Canopy\User\lib\site-packages\pandas\parser.pyd in pandas.parser.TextReader.__cinit__ (pandas\parser.c:3143)()

C:\Users\ngold\AppData\Local\Enthought\Canopy\User\lib\site-packages\pandas\parser.pyd in pandas.parser.TextReader._setup_parser_source (pandas\parser.c:5765)()

IOError: File 2014-10-20 0551_39_test1a_percol_rundata.csv does not exist



In [3]:

    
# 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 [4]:

    
# 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()









    



[1.0, 0.905, 0.819, 0.741, 0.67, 0.607, 0.549, 0.497, 0.449, 0.407]

Examining a Typical Run

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 [6]:

    
# Select the innovation sizes for one of the runs
innovSizes = invDF.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 = "TypicalRun_InnovHist.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 [7]:

    
# 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=(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 = "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 [8]:

    
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")









    



Optimization terminated successfully.
         Current function value: 0.377410
         Iterations 7
Optimization terminated successfully.
         Current function value: 0.377377
         Iterations 7
Optimization terminated successfully.
         Current function value: 0.377345
         Iterations 7
Optimization terminated successfully.
         Current function value: 0.377339
         Iterations 7
Optimization terminated successfully.
         Current function value: 0.247534
         Iterations 8
Optimization terminated successfully.
         Current function value: 0.247530
         Iterations 8
Optimization terminated successfully.
         Current function value: 0.247293
         Iterations 8
Optimization terminated successfully.
         Current function value: 0.247173
         Iterations 8

Simulation Tests - Plots and Analysis

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 [14]:

    
# 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 = "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 [15]:

    
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 = "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 [18]:

    
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 = "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 [19]:

    
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 = "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 [20]:

    
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 = "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 [21]:

    
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 = "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 [22]:

    
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 = "Test3B_MaxBPF.png".format(fileTime)
    plt.savefig(figName)
plt.show()

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 [3]:

    
fig = plt.figure(figsize=(13,5))
ax1 = fig.add_subplot(121)
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)

ax1 = fig.add_subplot(122)
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 = "{0}_microsim_maxbpf_resistmu.png".format(fileTime)
    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 [1]:

    
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.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()


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()









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-2445e0262351> in <module>()
      2 
      3 fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(3,4))
----> 4 dataPats = microSim[(microSim.doPatent==True)]['maxBPF']
      5 dataNoPats = microSim[(microSim.doPatent==False)]['maxBPF']
      6 data = np.array([dataPats,dataNoPats])

NameError: name 'microSim' is not defined



In [24]:

    
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.04)]['lmaxBPF'].plot(kind='kde',ax=axes,c='blue',lw=2)
patMean= mean(microSim[(microSim.doPatent==True)&(microSim.piDecayRate<0.04)]['lmaxBPF'])
microSim[(microSim.doPatent==False)&(microSim.piDecayRate<0.04)]['lmaxBPF'].plot(kind='kde',ax=axes,c='green',lw=2)
noPatMean= mean(microSim[(microSim.doPatent==False)&(microSim.piDecayRate<0.04)]['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.07)]['lmaxBPF'].plot(kind='kde',ax=axes,c='blue',lw=2)
patMean = mean(microSim[(microSim.doPatent==True)&(microSim.piDecayRate>0.07)]['lmaxBPF'])
microSim[(microSim.doPatent==False)&(microSim.piDecayRate>0.07)]['lmaxBPF'].plot(kind='kde',ax=axes,c='green',lw=2)
noPatMean = mean(microSim[(microSim.doPatent==False)&(microSim.piDecayRate>0.07)]['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)
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()