# pretty matplotlib figs
# load libraries and set plot parameters
import PrettyTable as pt
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'png')
plt.rcParams['savefig.dpi'] = 75

plt.rcParams['figure.autolayout'] = False
plt.rcParams['figure.figsize'] = 10, 6
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['axes.titlesize'] = 20
plt.rcParams['font.size'] = 16
plt.rcParams['lines.linewidth'] = 2.0
plt.rcParams['lines.markersize'] = 8
plt.rcParams['legend.fontsize'] = 14

plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.serif'] = "cm"
plt.rcParams['text.latex.preamble'] = r"\usepackage{subdepth}, \usepackage{type1cm}"

# pretty tables
data = np.array([[1,2,3],[2,3,4]])
pt.PrettyTable(data, [r"$\frac{a}{b}$", r"$b$", r"$c$"])
  # pt.PrettyTable takes a two-dimensional (!) numpy array and converts it to a table.

# easy reading of learge text files
def loadtxt(filename, ncols=None, dtyp=None, commentchars=["#","@"]):
    global data
    if ncols==None:
        with open(filename) as f:
            l = f.readline()
            while l[0] in commentchars:
                print(l)
                l = f.readline()
            ncols = len(l.split())
            data = np.fromfile(filename, sep=" ", dtype=dtyp)
    return data.reshape((len(data)/ncols, ncols))

# Difference betw .iloc & .loc (.ix has been deprecated):
#  Cf. http://pandas-docs.github.io/pandas-docs-travis/indexing.html#ix-indexer-is-deprecated
# .loc is primarily label based, but may also be used with a boolean array.
# .iloc is primarily integer position based (from 0 to length-1 of the axis), but may also be used with a boolean array.

# Below is an example
dfd = pd.DataFrame({'A': [1, 2, 3],\
                     'B': [4, 5, 6]},\
                    index=list('abc'))
# Deprecated way
# dfd.ix[[0, 2], 'A']

# Now recommended way
print(dfd.loc[dfd.index[[0, 2]], 'A'])  # this gives a pd.Series
print(dfd.loc[dfd.index[[0, 2]], ['A']])  # this gives a pd.DataFrame
print(dfd.iloc[[0, 2], dfd.columns.get_loc('A')])  # only applicable for ONE col; this gives a pd.Series
print(dfd.iloc[[0, 2], dfd.columns.get_indexer(['A', 'B'])])  # applicable for multicols; gives a pd.DataFrame

# Various max drawdown funcs:

# Func1: max drawdown for one col of ret
def mddoneret(ret):  
    # ret is a pd.Series
    r = ret.add(1).cumprod()
    dd = r.div(r.cummax()).sub(1)
    mdd = dd.min()
    end = dd.idxmin()
    start = r.loc[:end].idxmax()
    return mdd, start, end

# Func2: max active drawdown for two cols of ret(pf & bn)
def mddtworet(p, b):  
    p = p.add(1).cumprod()  # portfolio
    b = b.add(1).cumprod()  # benchmark
    pmb = p - b
    cam = pmb.expanding(min_periods=1).apply(lambda x: x.argmax())  # expanding window of cumulative argmax
    p0 = pd.Series(p.loc[cam.values.astype(int)].values, index=p.index)
    b0 = pd.Series(b.loc[cam.values.astype(int)].values, index=b.index)
    dd = p/p0 - b/b0
    mdd = dd.min()
    end = dd.idxmin()
    start = cam.loc[end]
    return mdd, start, end

# Func3: rolling max drawdown for one col of `ret` and for length of `win`
# Inputs: ret is pfret, win is window length (exogenous parameter)
# Outputs: tes is a pd.DataFrame, whose cols 'mdd, start, end' are 'maxdrawdown, mdd begday, mdd finday' resp
def rolmddoneret(ret, win):
    tes = pd.DataFrame(index=pd.date_range(start='20160101', periods=len(ret), freq='D'))
    tes['pfret'] = ret.values
    tes['pfnav'] = tes.pfret.add(1).cumprod()
    for i in range(len(tes)-win+1):
        rolp = tes.pfnav.iloc[i:i+win]
        cam = rolp.expanding(min_periods=1).apply(lambda x: x.argmax())  # expanding window
        rolp0 = pd.Series(rolp.iloc[cam.values.astype(int)].values, index=rolp.index)
        roldd = rolp/rolp0 - 1
        tes.loc[tes.index[i+win-1], 'mdd'] = roldd.min()
        tes.loc[tes.index[i+win-1], 'end'] = roldd.idxmin()
        tes.loc[tes.index[i+win-1], 'start'] = cam.index[int(cam.loc[tes.loc[tes.index[i+win-1], 'end']])]
    return tes

# Func4: rolling max active drawdown for two cols of ret (pf & bn)
# Inputs: p is pfret, b is bnret, (p,b must be of same length), win is window length (exogenous parameter)
# Outputs: tes is a pd.DataFrame, whose cols 'mdd, start, end' are same as mddtworet's meaning
def rolmddtworet(p, b, win):
    # tes = pd.DataFrame(index=pd.date_range(start='20160101', periods=200, freq='D'))
    # tes = tes.assign(pfret=p.values, bnret=b.values)  # `nav` is the new colname. Note: df.assign is NOT inplace!
    tes = pd.DataFrame(index=pd.date_range(start='20160101', periods=len(p), freq='D'))
    tes['pfret'] = p.values
    tes['bnret'] = b.values
    tes['pfnav'] = tes.pfret.add(1).cumprod()
    tes['bnnav'] = tes.bnret.add(1).cumprod()
    tes['pmb'] = tes.pfnav - tes.bnnav

    for i in range(len(tes)-win+1):
        rolp = tes.pfnav.iloc[i:i+win]
        rolb = tes.bnnav.iloc[i:i+win]
        rolpmb = tes.pmb.iloc[i:i+win]
        cam = rolpmb.expanding(min_periods=1).apply(lambda x: x.argmax())  # expanding window
        rolp0 = pd.Series(rolp.iloc[cam.values.astype(int)].values, index=rolp.index)
        rolb0 = pd.Series(rolb.iloc[cam.values.astype(int)].values, index=rolb.index)
        roldd = rolp/rolp0 - rolb/rolb0
        tes.loc[tes.index[i+win-1], 'mdd'] = roldd.min()
        tes.loc[tes.index[i+win-1], 'end'] = roldd.idxmin()
        tes.loc[tes.index[i+win-1], 'start'] = cam.index[int(cam.loc[tes.loc[tes.index[i+win-1], 'end']])]

    return tes

%matplotlib inline
np.random.seed(314)
p = pd.Series(np.random.randn(200) / 100 + 0.001)  # as pfret
b = pd.Series(np.random.randn(200) / 100 + 0.001)  # as bnret
win = 70  # window length

# tworets: p & b

mdd, sd, ed = mddtworet(p, b)
f, axarr = plt.subplots(2, 2, figsize=[16, 10], sharex='col', sharey='row')

tes = rolmddtworet(p, b, win)
sd2 = tes.index[int(sd)]
ed2 = tes.index[int(ed)]

tes[['pfnav', 'bnnav']].plot(title='Cumulative Net Asset Value (Portfolio vs Benchmark)', ax=axarr[0,0])
axarr[0,0].axvspan(sd2, ed2, alpha=0.1, color='r')  # Draw a vertical span (rectangle) from sd to ed.
axarr[0,0].legend(['Portfolio','Benchmark'])

tes[['pmb', 'mdd']].plot(title='Cumulative Active Return, Rolling \& Static Max Drawdown', ax=axarr[1,0])
axarr[1,0].axvspan(sd2, ed2, alpha=1, color='0.618')  # alpha(0-1) for transparency. Cf. matplotlib.patches.Polygon for kwargs
  # Gray shades can be given as `color=` a string encoding a float in the 0-1 range
axarr[1,0].axhline(y=mdd, color='r', linestyle='--', label='mdd')  # add a horizontal line of mdd
axarr[1,0].legend(['Active','Rolling MaxDD','Static MaxDD'])  # other loc=: best, upper/lower/center right/left/center, center

# oneret: p

mdd, sd, ed = mddoneret(p)
tes = rolmddoneret(p, win)    
sd2 = tes.index[int(sd)]
ed2 = tes.index[int(ed)]

tes[['pfnav']].plot(title='Cumulative Net Asset Value (Portfolio)', ax=axarr[0,1])
axarr[0,1].axvspan(sd2, ed2, alpha=0.1, color='r')  # Draw a vertical span (rectangle) from sd to ed.
axarr[0,1].legend(['Portfolio'])

tes[['mdd']].plot(title='Rolling \& Static Max Drawdown', ax=axarr[1,1])
axarr[1,1].axvspan(sd2, ed2, alpha=1, color='0.618')  # alpha(0-1) for transparency. Cf. matplotlib.patches.Polygon for kwargs
axarr[1,1].axhline(y=mdd, color='r', linestyle='--', label='mdd')  # add a horizontal line of mdd
axarr[1,1].legend(['Rolling MaxDD','Static MaxDD'], loc='lower right')  # other loc=: best, upper/lower/center right/left/center, center

# Fine-tune figure; show x ticks for top plots and y ticks for right plots
plt.setp([a.get_xticklabels() for a in axarr[0, :]], visible=True);  # semicolon to suppress meaningless output
plt.setp([a.get_yticklabels() for a in axarr[:, 1]], visible=True);

# # Below deprecated
# keys = ['Portfolio', 'Benchmark']
# cum = pd.concat([p, b], axis=1, keys=keys).add(1).cumprod()  # navs; keys in pd.concat is used as cols
# cum['Active'] = cum.Portfolio - cum.Benchmark  # exret
# cum[['Portfolio', 'Benchmark']].plot(title='Cumulative Absolute', ax=axarr[0,0])
# axarr[0,0].axvspan(sd, ed, alpha=0.1, color='r')  # Draw a vertical span (rectangle) from sd to ed. 
# cum[['Active']].plot(title='Cumulative Active', ax=axarr[1,0])
# axarr[1,0].axvspan(sd, ed, alpha=1, color='0.618')  # alpha(0-1) for transparency. Cf. matplotlib.patches.Polygon for kwargs
#   # Gray shades can be given as `color=` a string encoding a float in the 0-1 range

def previous_and_next(some_iterable):
    prevs, items, nexts = it.tee(some_iterable, 3)
    prevs = it.chain([None], prevs)
    nexts = it.chain(it.islice(nexts, 1, None), [None])
    return zip(prevs, items, nexts)
mylist = ['banana', 'orange', 'apple', 'kiwi', 'tomato']

for previous, item, nxt in previous_and_next(mylist):
    print ("Item is now", item, "next is", nxt, "previous is", previous)

# 求解线性方程组：x+3y+5z=10, 2x+5y+z=8, 2x+3y+8z=3;及将右侧三数换为(9,8,13)

A = np.mat('1,3,5; 2,5,1; 2,3,8')    # 构造系数矩阵 A
# B = np.mat('10,8,3; 9,8,13').T       # 构造矩阵(向量) b,注意尺寸(3*2)，否则不符合矩阵的乘法
# 上行与下行等价
B = np.array([[10,8,3],[9,8,13]]).T  # 构造二维数组 B
x = np.linalg.solve(A,B)             # 调用 solve 函数
print(x)  # x的每列是相应方程组的解。可用x.T[0],[1]的形式取出
print(type(x))

# 检验结果；下两行等价
print(A * x)
print(np.dot(A, x))

# 求解非线性方程组：x+y=4, x**2+y**2=8.
import scipy.optimize as so
def f(x):  # 预备求解形如f(x)=0的方程组
    x1, x2 = x[0], x[1]
    return [x1 + x2 - 4, x1**2 + x2**2 - 8]  # 返回误差
result = so.fsolve(f, [1,1])  # 第二参数为未知向量(x,y)的初始值
print ('the result is', result)
print ('the error is', f(result))

# 若方程组中未知数多，而各方程所含未知数少，Jacobi matrix稀疏，可传递Jacobi给fsolve以提高计算精度
def j(x):                         #定义雅可比矩阵
    x1,x2=x[0],x[1]
    return [[1,1],[2*x1,2*x2]]

result=so.fsolve(f,[1,1],fprime=j)   #带有雅可比矩阵的 fsolve 函数调用

print ('the result (Jacobi provided) is',result)
print ('the error (Jacobi provided) is',f(result))

# 推荐用此方法连接M$ SQL Server. Justyre 20180115
import pyodbc
import pandas as pd
# import getpass
# pwd = getpass.getpass(prompt='Enter the password: ')  # this is a password getter requiring manual input
pwd = '111111'
constr = pyodbc.connect(DRIVER='{SQL SERVER}', \
                        SERVER="192.168.1.126", DATABASE="WindDB", UID="wind", PWD=pwd, charset="UTF-8")
sqlcmd = '''select top 10 * from WindDB.dbo.ASHAREDESCRIPTION'''
mydata3 = pd.read_sql(sqlcmd, constr)
mydata3

data=pd.DataFrame({'price':[949.975 ,941.370,949.975,941.370,949.975],'volume':[0.01,0.01,0.01,0.01,0.01]})
data['timer']=['2014-01-15 14:29:54','2014-01-15 14:29:59','2014-01-15 14:30:17','2014-01-15 14:30:24','2014-01-15 14:30:36']
data.dtypes
data.timer=pd.to_datetime(data.timer)
data.set_index('timer',inplace=True)
data
a1 = data.price.resample('30s').ohlc()
a1['volsum'] = data.volume.resample('30s').sum()
a1

dat = pd.read_csv('Guerry.csv',index_col=0)
dat

# Fit regression model (using the natural log of one of the regressaors)
results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()

# Inspect the results
# print(results.summary())
results.summary()

# Explanation of output:

# TABLE 1
# Basics:
# No. Observations: Total number of observations used in regression, n=86.
# Df Residuals, Df Model: Total variance has df of n-1=85. 
#   Model has df of #(estimated coef)-1=3-1=2 (since intercept is included by default).
#   Thus, residual df = total df - model df = 85-2=83.

# Anova Table: analysis of variance.
#   Sum of Squares of Regression(ie. MOdel) (SSR) = \sum_{i=1}^n (Y_i^hat - Y_bar)^2（拟合值vs期望）; Mean SR=SSR/DfModel.
#   SS of Error(ie. Residuals) (SSE)              = \sum_{i=1}^n (Y_i - Y_i^hat)^2（原值vs拟合值）; Mean SE=SSE/DfResiduals.
#   SS of Total (SSTO)                            = \sum_{i=1}^n (Y_i - Y_bar)^2（原值vs期望）.

# Overall Model Fit:
# R-squared: Proportion of variance in `Lottery` which can be explained by the independent vars (ie. predictors).
# Adj R-squared: Penalizes addition of extraneous predictors to the model.
#   AdjRsq = 1 – ((1 – Rsq)((N – 1) /( N – k – 1)) where k is the number of predictors
# F-statistic: F(DfModel, DfResiduals) = MSR/MSE.
# Prob (F-statistic): Prob > F, used in testing the null hypothesis that all of the model coefficients are 0. Smaller, better.

# Model Quality for Model Selection:
# Log-Likelihood: The `likelihood` is the product of the density evaluated at the observations. Bigger, better.
  # Used to compare the fit of different coefs. Log-likelihood Lhat = Prob(x|thetahat, Model), thetahat is paraval that maxs L.
# AIC: Measure of relative quality of statmodels for a given set of data. Smaller, better.
  # AIC = 2 * #(estimated coef) - 2 * log(Lhat)
# BIC: Criterion for model selection. Smaller, better.
  # BIC = log(No. Observations) * #(estiamted coef) - 2 * log(Lhat)

# TABLE 2
# Parameter Estimates:
# coef: coefficients of predictors.
# stderr: standard errors associated with the coefficients.
# t: t-statistics used in testing whether a given coefficient is significantly different from zero.
# P>|t|: 2-tailed p-values used in testing the null hypothesis that the coefficient (parameter) is 0. Smaller, better.
# [0.025, 0.975]: the 95% confidence intervals for the coefficients.
#   The confidence intervals are related to the p-values such that 
#   the coefficient will not be statistically significant at alpha = .05 if the 95% confidence interval includes zero.

# TABLE 3  
# Omnibus: a set of tests like F-test
# Skewness: 0 for symmetric data, positive for right-skewed (right tail is longer than left tail).
# Kurtosis: 3 for std normdist, >3 for heavy-tailed distribution, <3 for light-tailed dist.
# Durbin-Watson: Detect autocorrelation in residuals. Always betw 0 and 4.
#   2 is no autocorr; <<2 (eg.<1) is evidence of positive serial corr; >>2 nega corr, ie. underesti of stat significance.
# Jarque-Bera: a goodness-of-fit test of whether sample data have the skewness and kurtosis matching a normal distribution.
#   Smaller, better. For normdist, JB=0.
# Condition number: measures how sensitive a function is to changes or errors in the output, 
#   and how much error in the output results from an error in the input. Lower, better.
#   Low: well-conditioned problem, High: ill-conditioned. Can be used as a diagnostic for multicollinearity.
#   If >30, regression has significant mulcoll.
#   Mulcoll exists if >=2 variables related to the high condition number have high proportions of variance explained.

dir(results)

# Normality test: Jarque-Bera
name = ['Jarque-Bera', 'Chi^2 two-tail prob.', 'Skew', 'Kurtosis']
test = sms.jarque_bera(results.resid)
list(zip(name, test))

# Omni test
name = ['Chi^2', 'Two-tail probability']
test = sms.omni_normtest(results.resid)
list(zip(name, test))

# Influence of each observation
from statsmodels.stats.outliers_influence import OLSInfluence
test_class = OLSInfluence(results)
test_class.dfbetas[:5,:]

dir(OLSInfluence)

%matplotlib inline
from statsmodels.graphics.regressionplots import plot_leverage_resid2
fig, ax = plt.subplots(figsize=(8,6))
fig = plot_leverage_resid2(results, ax = ax)

# Multicollinearity test
np.linalg.cond(results.model.exog), results.condition_number

# Heteroskedasticity test
# Breusch-Pagan test
name = ['Lagrange multiplier statistic', 'p-value',
        'f-value', 'f p-value']
test = sms.het_breuschpagan(results.resid, results.model.exog)
list(zip(name, test))

# Goldfeld-Quandt test
name = ['F statistic', 'p-value']
test = sms.het_goldfeldquandt(results.resid, results.model.exog)
list(zip(name, test))

# Linearity test
# Harvey-Collier multiplier test for Null hypothesis that the linear specification is correct
name = ['t value', 'p value']
test = sms.linear_harvey_collier(results)
list(zip(name, test))

results.params

all(results.resid == dat.Lottery - results.predict(dat))  # Validation: residual = observed - predicted

src = pd.concat([pd.Series([1] * len(dat), index=list(range(1, len(dat)+1))).rename('B0'), \
                 dat.Literacy, np.log(dat.Pop1831)], axis=1)  # prepares matrix X
src['Observed_y'] = dat.Lottery
src['Predicted_y'] = results.predict(dat)  # gives predicted values
src['Residual_y'] = results.resid  # residuals
src

# Lottery ~ Literacy + np.log(Pop1831)
plt.figure(figsize=(9, 4))
plt.subplot(121)
plt.scatter(src.Residual_y, src.Predicted_y, marker='o')
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Residual vs Predicted')
plt.subplot(122)
plt.scatter(src.Observed_y, src.Predicted_y, marker='o')
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Observed vs Predicted')

# A way to display Latex in cell output
from IPython.core.display import display, Latex
display(Latex(r'$\alpha_1$'))

import matplotlib.patches as mpat
import matplotlib.collections as mcol

fig, ax = plt.subplots(figsize=(3,3))
xs = [0, .5, .5,-.5,-.5]  # x of circle centers
ys = [0, .5,-.5, .5,-.5]  # y of circle centers
rs = [1,0.5,0.5,0.5,0.5]  # radius of circles
patches = []
for x, y, r in zip(xs, ys, rs):
    circle = mpat.Circle((x, y), r)
    patches.append(circle)

# patches.append(mpat.Rectangle((1,1),2,3))
p = mcol.PatchCollection(patches, alpha=0.4)
ax.add_collection(p)
ax.set_xlim(0,1)
ax.set_ylim(0,1)
ax.set_aspect('equal', adjustable='box')  # `ax` may be changed to `plt.gca()`

# Below deprecated
# fig1 = plt.figure()
# ax1 = fig1.add_subplot(111, aspect='equal')
# ax1.add_patch(mpat.Circle((0.5,0.5),1,color='b',fill=False))  # remove background

# circle1 = plt.Circle((0.5,0.5),0.5,color='r',fill=False)
# plt.gcf().gca().add_artist(circle1)  # gcf: get current figure, gca: get current axis

# fig1.savefig('rect1.png', dpi=90, bbox_inches='tight')

# 三弧问题的面积解法 by Justyre 20170415
# 首先，通过列两圆的参数方程，可解得交点的圆心角alpha
alpha = 1/2 * np.arcsin(3/4)
alpha, np.degrees(alpha)

# the two values of alpha are in [pi, 3/2 pi], [0, 1/2 pi], resp
a1 = 3/2 * np.pi - alpha
a2 = alpha
np.sin(2*a1), np.sin(2*a2), 3/4

theta = 1/2 * np.arcsin(-9/16)
theta, np.degrees(theta)

# the two values of theta are both in [3/2 pi, 2 pi]
t1 = 3/2 * np.pi - theta
t2 = 2 * np.pi + theta
np.sin(2*t1), np.sin(2*t2), -9/16

2 * np.cos(theta) - 1 - np.cos(alpha), 2 + 2 * np.sin(theta) - 1 - np.sin(alpha)

2 * np.cos(t1) - 1 - np.cos(a1), 2 + 2 * np.sin(t1) - 1 - np.sin(a1)

2 * np.cos(t2) - 1 - np.cos(a2), 2 + 2 * np.sin(t2) - 1 - np.sin(a2)

x1 = np.cos(t1)
y1 = 1 + np.sin(t1)
x2 = np.cos(t2)
y2 = 1 + np.sin(t2)

x1p = 1/2 * (1 + np.cos(a1))  # should be same as x1
y1p = 1/2 * (1 + np.sin(a1))
x2p = 1/2 * (1 + np.cos(a2))
y2p = 1/2 * (1 + np.sin(a2))

x1, x1p, y1, y1p, x2, x2p, y2, y2p

# 扇形面积
angle = 2 * a2 + np.pi / 2  # 圆心角
shan = angle / (2 * np.pi) * np.pi * (1/2 ** 2)  # 扇形面积
shan

import scipy.integrate as si
# 左边方程
s1 = (y1 - 1/2) / (x1 - 1/2)  # slope
i1 = 1/2 - 1/2 * s1  # intercept
# 右边方程
s2 = (y2 - 1/2) / (x2 - 1/2)  # slope
i2 = 1/2 - 1/2 * s2  # intercept
print(s1, i1, s2, i2)

# 左边下面积：int_x1^0.5 (s1 x + i1) dx
l1area = si.quad(lambda x: s1 * x + i1, x1, 1/2)[0]
# 左边下圆弧下面积
c1area = si.quad(lambda x: 1 - np.sqrt(1-x**2), x1, 1/2)[0]
# 右边下面积
l2area = si.quad(lambda x: s2 * x + i2, 1/2, x2)[0]
# 右边下圆弧下面积
c2area = si.quad(lambda x: 1 - np.sqrt(1-x**2), 1/2, x2)[0]
print(l1area, c1area, l2area, c2area)

# 弧边三角形面积
curtriangle = l1area - c1area + l2area - c2area
curtriangle

# 所求阴影面积
(shan - curtriangle) * 2

def f(x):
    return x - x/2*np.sqrt(1-x**2) - 1/2*np.arcsin(x)
print(f(1/2) - f(x1), f(x2) - f(1/2))  # 验证此为上述c1,c2积分的有限形式
def f(x):
    return s1*x**2/2 + i1*x
print(f(1/2) - f(x1))  # 验证此为上述l1积分的有限形式

# 五类账户的方向：资产+成本=负债+权益+损益（损左益右），AC=DER

# 例2-1：追溯调整法
a = 450e4  # 05购入成本
b = 110e4  # 06购入成本
t = .25  # 所得税率
法定盈余公积比例 = .1
任意盈余公积比例 = .05
普通股数 = 4500e4
a05 = 510e4  # a05年末公允价值
a06 = 510e4
b06 = 130e4
税前差异05 = a05 - a  # 05期末06期初
累积影响数05 = 税前差异05 * (1 - t)  # 05期末06期初
税前差异06 = a06 + b06 - a05 - b
累积影响数06 = 税前差异06 * (1 - t)
税前差异05, 累积影响数05, 税前差异06, 累积影响数06

# 对05年事项调整分录

# 调整累积影响数
print('20051231')
print('借：交易性金融资产——公允价值变动\t', 税前差异05)  # 资产+
print(q+'贷：利润分配——未分配利润\t', 累积影响数05)  # 益+
print(q+q+'递延所得税负债\t\t', 税前差异05 * t)  # 负债+
# 调整利润分配
print('借：利润分配——未分配利润\t\t', round(累积影响数05 * (法定盈余公积比例 + 任意盈余公积比例),2))  # 益-
print(q+'贷：盈余公积——法定盈余公积\t', 累积影响数05 * 法定盈余公积比例)  # 益+
print(q+q+'盈余公积——任意盈余公积\t', 累积影响数05 * 任意盈余公积比例)  # 益+

# 06年事项

# 调整累积影响数
print('20061231')
print('借：交易性金融资产——公允价值变动\t', 税前差异06)  # 资产+
print(q+'贷：利润分配——未分配利润\t', 累积影响数06)  # 益+
print(q+q+'递延所得税负债\t\t', 税前差异06 * t)  # 负债+
# 调整利润分配
print('借：利润分配——未分配利润\t\t', round(累积影响数06 * (法定盈余公积比例 + 任意盈余公积比例),2))  # 益-
print(q+'贷：盈余公积——法定盈余公积\t', 累积影响数06 * 法定盈余公积比例)  # 益+
print(q+q+'盈余公积——任意盈余公积\t', 累积影响数06 * 任意盈余公积比例)  # 益+

# 财务报表重述。站在07年初

print('资产负债表的调整')
delta = {'以公允价值计量且其变动计入当期损益的金融资产07年初余额+': 600000.0 + 200000.0, \
         '递延所得税负债+': 150000.0 + 50000.0, \
         '盈余公积+': 67500.0 + 22500.0, \
         '利润分配——未分配利润（即留存收益）07年初余额+': 450000.0 - 67500.0 + 150000.0 - 22500.0
        }
print(delta)

print('利润表的调整')
delta = {'公允价值变动收益06年金额+': 200000.0, \
         '所得税费用06年金额+': 50000.0, \
         '净利润06年金额+': 150000.0, \
         '基本每股收益06年金额+': round(150000.0 / 普通股数, 4)
        }
print(delta)

print('所有者权益变动表的调整')
delta = {'盈余公积06年初金额+': 67500.0, \
         '未分配利润06年初金额+': 450000.0 - 67500.0, \
         '所有者权益合计06年初金额+': 450000.0, \
         '盈余公积06年金额+': 22500.0, \
         '未分配利润06年金额+': 150000.0 - 22500.0, \
         '所有者权益合计06年金额+': 150000.0, \
         '盈余公积07年初金额+': 67500.0 + 22500.0, \
         '未分配利润07年初金额+': 450000.0 - 67500.0 + 150000.0 - 22500.0, \
         '所有者权益合计07年初金额+': 450000.0 + 150000.0
        }
print(delta)

# 例2-2：未来适用法
期初存货 = 250e4
购入存货实际成本 = 1800e4
期末存货先进先出法 = 450e4
当年销售额 = 2500e4
当年其他费用 = 120e4
t = .25
期末存货后进先出法 = 220e4

print('后进先出法营业成本', 期初存货 + 购入存货实际成本 - 期末存货后进先出法)
print('先进先出法营业成本', 期初存货 + 购入存货实际成本 - 期末存货先进先出法)
print('改用FIFO后，净利润增加', (期末存货先进先出法 - 期末存货后进先出法) * (1-t))

# 存货分类：原材料、在产品、半成品、产成品、商品、周转材料
# 存货成本 = 采购 + 加工 + 使存货达到目前场所和状态所发生的其他成本
# 外购存货的成本 = 采购成本 = 购买价款 + 税费 + 运输装卸保险等费用
#   已销售商品的进货费用 → 主营业务成本；未售商品的进货费用 → 期末存货成本；较小的 → 销售费用

# 存货周转的相关分录：
#   借：生产成本（直接构成产品实体的成本）、制造费用（成本类，间接费用）、管理费用、销售费用，贷：原材料、应付职工薪酬、累计折旧、银行存款。
#   借：生产成本，贷：制造费用。（表示未完工产品成本）
#   借：库存商品，贷：生产成本。（表示已完工产品成本）
#   销售时：借：主营业务成本，贷：库存商品。
# 投资者投入存货的成本：协议价优先，不公允的按公允价值。
#   分录：借：库存商品、应交税费——应交增值税（进项税额），贷：实收资本、资本公积——资本溢价。
# 盘盈存货的成本：按重置成本，经待处理财产损溢，管理报批后冲减当期管理费用
# 采购用于广告营销的商品、向客户预付款未取得商品时：预付账款，取得商品后入销售费用

# 发出存货的计量：FIFO，移动加权平均，月末一次加权平均，个别计价。
# 销售存货的分录：借：营业成本（主营业务成本or其他业务成本） + 存货跌价准备，贷：库存商品or原材料。
# 企业周转材料的领用分录：借：成本费用，贷：周转材料——包装物or低值易耗品。

# 待处理财产损溢（损益类）：借方为盘亏，贷方为盘盈。期末处理后，该科目无余额
# 存货盘盈的分录：
#   借：原材料or库存商品，贷：待处理财产损溢。
#   报批后：借：待处理财产损溢，贷：管理费用。
# 存货盘亏的分录：
#   借：待处理财产损溢，贷：原材料or库存商品 + （购进存货特有）应交税金──应交增值税(进项税额转出)。
#   借：管理费用（源于计量收发差错和管理不善等）or其他应收款（源于过失人赔偿、保险赔偿）or营业外支出（源于自然灾害），贷：待处理财产损溢。

# 例4-1
# 购入需要安装的固定资产

print('20090201')
# 支付设备价款、增值税、运输费
print('借：在建工程——某设备\t\t', 50e4 + 2500 * 0.93)  # 运输费应计入在建工程成本。deprecated: 运输费价内增值税（现行一般都是价外了）
print(q+'应交税费——应交增值税（进项税额）', 8.5e4 + 2500 * 0.07)  # 假定运输费增值税率7%
print(q+'贷：银行存款\t\t\t', 50e4 + 8.5e4 + 2500)
# 领用公司原材料、支付安装工人工资
print('借：在建工程——某设备\t\t', 3e4 + 4900)
print(q+'贷：原材料\t\t\t', 3e4)  # 购进原材料时的增值税在“售出”原材料时才抵扣，“使用”原材料时用不到
print(q+q+'应付职工薪酬\t\t', 4900.0)
# 设备安装完毕，到达预定可使用状态
print('借：固定资产——某设备\t\t', 50e4 + 2500 * 0.93 + 3e4 + 4900)
print(q+'贷：在建工程——某设备\t\t', 50e4 + 2500 * 0.93 + 3e4 + 4900)

# 例4-2
# 超出正常信用条件下购买固定资产，具有融资性质。涉及未确认融资费用的分摊
总价 = 900e4  # 不考虑增值税
期数 = 5 * 2
折现率 = 0.1  # 每期折现率

年金 = 总价 / 期数  # 每期付款
折现系数PA = round(np.pv(rate=折现率, nper=期数, pmt=-1, fv=0, when='end'), 4)  # numpy.pv(rate, nper, pmt, fv=0.0, when='end')
现值 = round(年金 * 折现系数PA, 2)
年金, 折现系数PA, 现值

# 购买设备

print('20170101')
print('借：在建工程——某设备\t', 现值)  # 总本金（资产+）
在建工程 = 现值  # 供结转用
print(q+'未确认融资费用\t', 总价 - 现值)  # 总利息（损+）
print(q+'贷：长期应付款——乙公司\t', 总价)  # 总本息和（负债+）

# 20170101-20171231是设备安装期间，未确认融资费用的分摊额可计入固定资产成本

print('20170630')
剩余本金 = 现值  # 本次归还之前的剩余本金
本次融资费用 = round(剩余本金 * 折现率, 2)  # 未确认融资费用（即利息）在本次变成已确认的部分
print('借：在建工程——某设备\t', 本次融资费用)  # 资产+
在建工程 += 本次融资费用  # 供结转用
print(q+'贷：未确认融资费用\t', 本次融资费用)  # 损-
print('借：长期应付款——乙公司\t', 年金)  # 负债-：年金同时归还了本金和利息
print(q+'贷：银行存款\t\t', 年金)  # 资产-

print('20171231')
剩余本金 -= 年金 - 本次融资费用  # 本次归还之前的剩余本金，这里的“本次融资费用”是上次发生的
本次融资费用 = round(剩余本金 * 折现率, 2)  # 未确认融资费用（即利息）在本次变成已确认的部分
print('借：在建工程——某设备\t', 本次融资费用)  # 资产+
在建工程 += 本次融资费用  # 供结转用
print(q+'贷：未确认融资费用\t', 本次融资费用)  # 损-
print('借：长期应付款——乙公司\t', 年金)  # 负债-：年金同时归还了本金和利息
print(q+'贷：银行存款\t\t', 年金)  # 资产-
# 本次特有分录：安装费，计入固定资产成本
安装费 = 398530.60
print('借：在建工程——某设备\t', 安装费)
在建工程 += 安装费  # 供结转用
print(q+'贷：银行存款\t\t', 安装费)
# 本次特有分录：已达到预定可使用状态，结转在建工程
print('借：固定资产——某设备\t', 在建工程)  # 资产+
print(q+'贷：在建工程——某设备\t', 在建工程)  # 资产-

# 20180101-20211231期间，该设备已达到预定可使用状态，未确认融资费用的分摊额应计入当期损益（而不再符合资本化条件）

ys = range(2018,2022)  # years: range(start,end,step)
ds = ['0630','1231']  # dates
for y, d in it.product(ys, ds):  # flat is better than nested
    print(str(y)+d)
    剩余本金 -= 年金 - 本次融资费用  # 本次归还之前的剩余本金，这里的“本次融资费用”是上次发生的
    print(剩余本金)
    本次融资费用 = round(剩余本金 * 折现率, 2)  # 未确认融资费用（即利息）在本次变成已确认的部分
    if str(y)+d == '20211231':  # 由于折现系数PA只保留了4位小数，末期的融资费用需要进行尾数调整
        本次融资费用 = 年金 - 剩余本金  # 年金同时归还了本金和利息，末期应将本金结清，所余部分即是剩余利息 
    print('借：财务费用\t\t', 本次融资费用)  # 损+
    print(q+'贷：未确认融资费用\t', 本次融资费用)  # 损-
    print('借：长期应付款——乙公司\t', 年金)  # 负债-
    print(q+'贷：银行存款\t\t', 年金)  # 资产-

# 高危行业提取安全生产费的分录：借：相关产品成本，贷：专项储备——安全生产费。
# 高危行业，使用提取的安全生产费形成固定资产的分录：（专项储备为权益类）
#   借：在建工程，贷：各项支出。（达到预定可使用状态时）借：固定资产a，贷：在建工程a。借：专项储备——安全生产费a，贷：累计折旧a。

# 出包方式建造固定资产的成本：建筑工程支出、安装工程支出、待摊支出。
# 分录：
#  施工期间：借：在建工程——建筑工程——某工程、在建工程——安装工程——某工程，贷：银行存款、预付账款
#  工程完成时，补付工程款：借：在建工程，贷：银行存款
#  将需安装设备运抵现场安装时：借：在建工程——在安装设备——某设备，贷：工程物资——某设备
#  为建造固定资产发生的待摊支出：借：在建工程——待摊支出，贷：银行存款、应付职工薪酬
# 在建工程达到预定可使用状态时，先将待摊支出按各工程的(建筑工程+安装工程+在安装设备)进行分配；然后如下分录：
#  借：固定资产，贷：在建工程——（建筑工程or安装工程or待摊支出）

# 盘盈固定资产：视为前期差错，经以前年度损益调整核算

# 例4-3
# 存在弃置费用的固定资产（仅限特殊行业的特定固定资产）
建造成本 = 250e8  # 元
预计使用寿命 = 40  # 年
弃置费用 = 25e8  # 将该值折现至今，即为预计负债
折现率 = 0.1
PF = np.pv(rate=折现率, nper=预计使用寿命, pmt=0, fv=-1, when='end')  # numpy.pv(rate, nper, pmt, fv=0.0, when='end')
PF

# 启用时
print('20170101')
print('借：固定资产——某核反应堆\t\t\t', 建造成本 + 弃置费用 * PF)
print(q+'贷：在建工程——某核反应堆\t\t', 建造成本)
print(q+q+'预计负债——某核反应堆——弃置费用\t', 弃置费用 * PF)
剩余本金 = 弃置费用 * PF

# 此后每一年
ys = range(2017,2057)  # years: range(start,end,step)
ds = ['1231']  # dates
for y, d in it.product(ys, ds):  # flat is better than nested
    print(str(y)+d)
    本次财务费用 = 剩余本金 * 折现率
    print('借：财务费用\t\t\t\t', 本次财务费用)
    print(q+'贷：预计负债——某核反应堆——弃置费用\t', 本次财务费用)
    print('已完成预计负债：', 剩余本金 + 本次财务费用)  # 即从启用时至今贷方预计负债的累计值
    剩余本金 *= 1 + 折现率

# 固定资产折旧：(]。当月增加的当月不计提，当月减少的当月仍计提。比较：无形资产摊销：[)。【折开摊闭】
# 双倍余额递减法不考虑预计净残值，因而在折旧到期前2年时要改为直线法（考虑净残值了）。
# 例4-4：双倍余额递减法
原价 = 120
预计使用寿命 = 5
预计净残值率 = 0.04
应提折旧01 = 原价 * 2 / 预计使用寿命  # 第一年
应提折旧02 = (原价 - 应提折旧01) * 2 / 预计使用寿命
应提折旧03 = (原价 - 应提折旧01 - 应提折旧02) * 2 / 预计使用寿命
# 末两年改用直线法
应提折旧04 = (原价 * (1 - 预计净残值率) - 应提折旧01 - 应提折旧02 - 应提折旧03) / 2
应提折旧05 = 应提折旧04
应提折旧01, 应提折旧02, 应提折旧03, 应提折旧04, 应提折旧05