This code snippet is a part of stackexchange answer http://bitcoin.stackexchange.com/a/20060/2343

solo.py

#!/usr/bin/env python
#-*- encoding: utf-8 -*-

# This code is a part of stackexchange answer
# http://bitcoin.stackexchange.com/a/20060/2343
# License: cc-wiki with attribution required

import numpy as np

HR = 550 # GHash/s, own hashrate
THR0 = 1.2e7 # GHash/s, network hashrate
DA = 1.2 # Expected difficulty adjustment
N = 26*2016 # For how much block make prediction? ~ 1 year

NB = 2016 # Number of block between difficulty adjustments
T0 = 10 # minutes, desired block time

def THR(n):
    return THR0*np.exp((n/NB)*np.log(DA))

p = HR/THR(np.arange(N))

# To caclulate distribution we use algorithm from
# Yili Hong "On Computing the Distribution Function for the Sum of
# Independent and Non-identical Random Indicators", 2011
# http://www.stat.vt.edu/research/Technical_Reports/TechReport11-2.pdf

x = np.zeros(N+1, dtype=np.complex128)
w = 2*np.pi / (N+1)
x[0] = 1.0 + 0.j
for l in xrange(1, int(np.ceil(N/2))+1):
    z = 1.0 - p[:] + p[:]*np.cos(w*l) + 1.j * p[:]*np.sin(w*l)
    dl = np.exp(np.sum(np.log(np.absolute(z[:]))))
    ang = np.sum(np.angle(z[:]))
    al = dl*np.cos(ang)
    bl = dl*np.sin(ang)
    x[l]     = al + 1.j*bl
    if l%1000 == 0: print l, '/', N/2

for l in xrange(int(np.ceil(N/2))+1, N+1):
    x[l] = np.conj(x[N+1-l])

x = x[:] / (N+1)
y = np.fft.fft(x)

print DA
for i in enumerate(np.real(y[:10])):
    print i[0], i[1]