Solve Wooldridge exercises in Python

Introduction

--Recently, I started studying econometrics, so I will solve the exercises in the textbook with Python.

• https://github.com/lystahi/wooldridge_python_notebook ――Since it is still in the middle, I think there are some parts that have not been completed or are wrong.

About textbooks

• Introductory Econometrics: A Modern Approach, by Jeffrey M. Wooldridge.
• https://www.amazon.co.jp/dp/1337558869/ --The solution by R can be found at here.

Solve

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
%matplotlib inline

import warnings
warnings.simplefilter('ignore')

Introduction

This notebook contains examples from Introductory Econometrics: A Modern Approach, 6e by Jeffrey M. Wooldridge. Each example illustrates how to load data, build econometric models, and compute estimates with Python.

Now, install and load the wooldridge package and lets get started!

# !pip install wooldridge

from wooldridge import *

Chapter 2: The Simple Regression Model

Example 2.10: A Log Wage Equation

Load the wage1 data and check out the documentation.

dataWoo('wage1', description=True)

name of dataset: wage1
no of variables: 24
no of observations: 526

+----------+---------------------------------+
 | variable | label                           |
+----------+---------------------------------+
 | wage     | average hourly earnings         |
 | educ     | years of education              |
 | exper    | years potential experience      |
 | tenure   | years with current employer     |
 | nonwhite | =1 if nonwhite                  |
 | female   | =1 if female                    |
 | married  | =1 if married                   |
 | numdep   | number of dependents            |
 | smsa     | =1 if live in SMSA              |
 | northcen | =1 if live in north central U.S |
 | south    | =1 if live in southern region   |
 | west     | =1 if live in western region    |
 | construc | =1 if work in construc. indus.  |
 | ndurman  | =1 if in nondur. manuf. indus.  |
 | trcommpu | =1 if in trans, commun, pub ut  |
 | trade    | =1 if in wholesale or retail    |
 | services | =1 if in services indus.        |
 | profserv | =1 if in prof. serv. indus.     |
 | profocc  | =1 if in profess. occupation    |
 | clerocc  | =1 if in clerical occupation    |
 | servocc  | =1 if in service occupation     |
 | lwage    | log(wage)                       |
 | expersq  | exper^2                         |
 | tenursq  | tenure^2                        |
+----------+---------------------------------+

These are data from the 1976 Current Population Survey, collected by
Henry Farber when he and I were colleagues at MIT in 1988.

The documentation indicates these are data from the 1976 Current Population Survey, collected by Henry Farber when he and Wooldridge were colleagues at MIT in 1988.

wage: average hourly earnings

educ: years of education

First, make a scatter-plot of the two variables and look for possible patterns in the relationship between them.

wage1 = dataWoo('wage1')

plt.scatter(wage1.educ, wage1.wage)
plt.title("Wages vs. Education, 1976")
plt.xlabel("years of education")
plt.ylabel("Hourly wages")
plt.show()

It appears that on average, more years of education, leads to higher wages.

The example in the text investigates what the percentage change between wages and education might be. So, we must use the log(wage).

Build a linear model to estimate the relationship between the log of wage (lwage) and education (educ).

y = wage1.lwage
x = wage1.educ

model = sm.OLS(x, y)
results = model.fit()

Print the summary of the results.

print(results.summary())

                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:                   educ   R-squared (uncentered):                   0.916
Model:                            OLS   Adj. R-squared (uncentered):              0.916
Method:                 Least Squares   F-statistic:                              5717.
Date:                Sat, 09 May 2020   Prob (F-statistic):                   2.18e-284
Time:                        15:36:54   Log-Likelihood:                         -1438.9
No. Observations:                 526   AIC:                                      2880.
Df Residuals:                     525   BIC:                                      2884.
Df Model:                           1                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                  coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
lwage          7.2081      0.095     75.608      0.000       7.021       7.395
==============================================================================
Omnibus:                       12.119   Durbin-Watson:                   1.699
Prob(Omnibus):                  0.002   Jarque-Bera (JB):               13.454
Skew:                          -0.295   Prob(JB):                      0.00120
Kurtosis:                       3.515   Cond. No.                         1.00
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Plot the log(wage) vs educ, adding a line representing the least squares fit.

x_with_const = sm.add_constant(x)

model = sm.OLS(x_with_const, y)
results = model.fit()

slope, intercept = results.params[1].values/100, results.params[0].values


def abline(slope, intercept):
    """Plot a line from slope and intercept"""
    axes = plt.gca()
    x_vals = np.array(axes.get_xlim())
    y_vals = intercept + slope * x_vals
    plt.plot(x_vals, y_vals, '--')


plt.scatter(x, y)
abline(slope, intercept)
plt.title("A Log Wage Equation")
plt.xlabel("years of education")
plt.ylabel("log of average hourly wages")
plt.show()

Chapter 3: Multiple Regression Analysis: Estimation

Example 3.2: Hourly Wage Equation

Check the documentation for variable information

dataWoo('wage1', description=True)

name of dataset: wage1
no of variables: 24
no of observations: 526

+----------+---------------------------------+
 | variable | label                           |
+----------+---------------------------------+
 | wage     | average hourly earnings         |
 | educ     | years of education              |
 | exper    | years potential experience      |
 | tenure   | years with current employer     |
 | nonwhite | =1 if nonwhite                  |
 | female   | =1 if female                    |
 | married  | =1 if married                   |
 | numdep   | number of dependents            |
 | smsa     | =1 if live in SMSA              |
 | northcen | =1 if live in north central U.S |
 | south    | =1 if live in southern region   |
 | west     | =1 if live in western region    |
 | construc | =1 if work in construc. indus.  |
 | ndurman  | =1 if in nondur. manuf. indus.  |
 | trcommpu | =1 if in trans, commun, pub ut  |
 | trade    | =1 if in wholesale or retail    |
 | services | =1 if in services indus.        |
 | profserv | =1 if in prof. serv. indus.     |
 | profocc  | =1 if in profess. occupation    |
 | clerocc  | =1 if in clerical occupation    |
 | servocc  | =1 if in service occupation     |
 | lwage    | log(wage)                       |
 | expersq  | exper^2                         |
 | tenursq  | tenure^2                        |
+----------+---------------------------------+

These are data from the 1976 Current Population Survey, collected by
Henry Farber when he and I were colleagues at MIT in 1988.

lwage: log of the average hourly earnings

educ: years of education

exper: years of potential experience

tenutre: years with current employer

Plot the variables against lwage and compare their distributions and slope (\beta) of the simple regression lines.

y = wage1.lwage
x = wage1.educ

model = sm.OLS(x, y)
results = model.fit()

x_with_const = sm.add_constant(x)

model = sm.OLS(x_with_const, y)
results = model.fit()

slope, intercept = results.params[1].values/100, results.params[0].values

plt.scatter(x, y)
abline(slope, intercept)
plt.title("A Log Wage Equation")
plt.xlabel("years of education")
plt.ylabel("log of average hourly wages")
plt.show()

y = wage1.lwage
x = wage1.exper

model = sm.OLS(x, y)
results = model.fit()

x_with_const = sm.add_constant(x)

model = sm.OLS(x_with_const, y)
results = model.fit()

slope, intercept = results.params[1].values/100, results.params[0].values

plt.scatter(x, y)
abline(slope, intercept)
plt.title("A Log Wage Equation")
plt.xlabel("years of education")
plt.ylabel("log of average hourly wages")
plt.show()

y = wage1.lwage
x = wage1.tenure

model = sm.OLS(x, y)
results = model.fit()

x_with_const = sm.add_constant(x)

model = sm.OLS(x_with_const, y)
results = model.fit()

slope, intercept = results.params[1].values/100, results.params[0].values

plt.scatter(x, y)
abline(slope, intercept)
plt.title("A Log Wage Equation")
plt.xlabel("years of education")
plt.ylabel("log of average hourly wages")
plt.show()

Estimate the model regressing educ, exper, and tenure against log(wage).

results = smf.ols("lwage ~ educ + exper + tenure", data=wage1).fit()
results.params

Intercept    0.284360
educ         0.092029
exper        0.004121
tenure       0.022067
dtype: float64

coef_df = pd.DataFrame(results.params)
coef_df.drop(index=['Intercept'], inplace=True)
coef_df.rename(columns={0: "coef"}, inplace=True)

coef_df

coef_df.sort_values(by='coef', ascending=False).plot.bar()

<matplotlib.axes._subplots.AxesSubplot at 0x1c1b4c9390>

Chapter 4: Multiple Regression Analysis: Inference

Example 4.1 Hourly Wage Equation

Using the same model estimated in example: 3.2, examine and compare the standard errors associated with each coefficient. Like the textbook, these are contained in parenthesis next to each associated coefficient.

results = smf.ols("lwage ~ educ + exper + tenure", data=wage1).fit()
results.summary()

results.params[2]/results.bse[2]

2.3914371073142298

Example 4.7 Effect of Job Training on Firm Scrap Rates

Load the jtrain data set.

dataWoo('jtrain', description=True)

name of dataset: jtrain
no of variables: 30
no of observations: 471

+----------+---------------------------------+
 | variable | label                           |
+----------+---------------------------------+
 | year     | 1987, 1988, or 1989             |
 | fcode    | firm code number                |
 | employ   | # employees at plant            |
 | sales    | annual sales, $                 |
 | avgsal   | average employee salary         |
 | scrap    | scrap rate (per 100 items)      |
 | rework   | rework rate (per 100 items)     |
 | tothrs   | total hours training            |
 | union    | =1 if unionized                 |
 | grant    | = 1 if received grant           |
 | d89      | = 1 if year = 1989              |
 | d88      | = 1 if year = 1988              |
 | totrain  | total employees trained         |
 | hrsemp   | tothrs/totrain                  |
 | lscrap   | log(scrap)                      |
 | lemploy  | log(employ)                     |
 | lsales   | log(sales)                      |
 | lrework  | log(rework)                     |
 | lhrsemp  | log(1 + hrsemp)                 |
 | lscrap_1 | lagged lscrap; missing 1987     |
 | grant_1  | lagged grant; assumed 0 in 1987 |
 | clscrap  | lscrap - lscrap_1; year > 1987  |
 | cgrant   | grant - grant_1                 |
 | clemploy | lemploy - lemploy[_n-1]         |
 | clsales  | lavgsal - lavgsal[_n-1]         |
 | lavgsal  | log(avgsal)                     |
 | clavgsal | lavgsal - lavgsal[_n-1]         |
 | cgrant_1 | cgrant[_n-1]                    |
 | chrsemp  | hrsemp - hrsemp[_n-1]           |
 | clhrsemp | lhrsemp - lhrsemp[_n-1]         |
+----------+---------------------------------+

H. Holzer, R. Block, M. Cheatham, and J. Knott (1993), “Are Training
Subsidies Effective? The Michigan Experience,” Industrial and Labor
Relations Review 46, 625-636. The authors kindly provided the data.

jtrain = dataWoo('jtrain')
jtrain_subset = jtrain[jtrain.year == 1987][['year', 'union', 'lscrap',
                                             'hrsemp', 'lsales', 'lemploy']]
jtrain_subset.head()

jtrain_subset.isnull().sum()

year         0
union        0
lscrap     103
hrsemp      28
lsales      38
lemploy     13
dtype: int64

jtrain_clean = jtrain_subset.dropna()

Now create the linear model regressing hrsemp(total hours training/total employees trained), lsales(log of annual sales), and lemploy(the log of the number of the employees), against lscrap(the log of the scrape rate).

results = smf.ols("lscrap ~ hrsemp + lsales + lemploy",
                  data=jtrain_clean).fit()
results.summary()

Chapter 5: Multiple Regression Analysis: OLS Asymptotics

Example 5.1: Housing Prices and Distance From an Incinerator

Load the hprice3 data set.

dataWoo('hprice3', description=True)

name of dataset: hprice3
no of variables: 19
no of observations: 321

+----------+---------------------------------+
 | variable | label                           |
+----------+---------------------------------+
 | year     | 1978, 1981                      |
 | age      | age of house                    |
 | agesq    | age^2                           |
 | nbh      | neighborhood, 1-6               |
 | cbd      | dist. to cent. bus. dstrct, ft. |
 | inst     | dist. to interstate, ft.        |
 | linst    | log(inst)                       |
 | price    | selling price                   |
 | rooms    | # rooms in house                |
 | area     | square footage of house         |
 | land     | square footage lot              |
 | baths    | # bathrooms                     |
 | dist     | dist. from house to incin., ft. |
 | ldist    | log(dist)                       |
 | lprice   | log(price)                      |
 | y81      | =1 if year = 1981               |
 | larea    | log(area)                       |
 | lland    | log(land)                       |
 | linstsq  | linst^2                         |
+----------+---------------------------------+

-

hprice3 = dataWoo('hprice3')

Next, model the log(price) against the log(dist) to estimate the percentage relationship between the two.

price_dist_model = smf.ols("lprice ~ ldist", data=hprice3).fit()

Create another model that controls for "quality" variables, such as square footage area per house.

price_area_model = smf.ols("lprice ~ ldist + larea", data=hprice3).fit()

price_dist_model.summary()

price_area_model.summary()

Chapter 6: Multiple Regression: Further Issues

Example 6.1: Effects of Pollution on Housing Prices, standardized.

Load the hprice2 data and view the documentation.

dataWoo('hprice2', description=True)

name of dataset: hprice2
no of variables: 12
no of observations: 506

+----------+-------------------------------+
 | variable | label                         |
+----------+-------------------------------+
 | price    | median housing price, $       |
 | crime    | crimes committed per capita   |
 | nox      | nit ox concen; parts per 100m |
 | rooms    | avg number of rooms           |
 | dist     | wght dist to 5 employ centers |
 | radial   | access. index to rad. hghwys  |
 | proptax  | property tax per $1000        |
 | stratio  | average student-teacher ratio |
 | lowstat  | perc of people 'lower status' |
 | lprice   | log(price)                    |
 | lnox     | log(nox)                      |
 | lproptax | log(proptax)                  |
+----------+-------------------------------+

D. Harrison and D.L. Rubinfeld (1978), “Hedonic Housing Prices and the
Demand for Clean Air,” by Harrison, D. and D.L.Rubinfeld, Journal of
Environmental Economics and Management 5, 81-102. Diego Garcia, a
former Ph.D. student in economics at MIT, kindly provided these data,
which he obtained from the book Regression Diagnostics: Identifying
Influential Data and Sources of Collinearity, by D.A. Belsey, E. Kuh,
and R. Welsch, 1990. New York: Wiley.

Estimate the usual lm model.

hprice2 = dataWoo('hprice2')
housing_level = smf.ols("price ~ nox + crime + rooms + dist + stratio",
                        data=hprice2).fit()

Estimate the same model, but standardized coefficients by wrapping each variable with R’s scale function:

from sklearn.preprocessing import scale

zprice = pd.DataFrame({"zprice":scale(hprice2["price"])})
znox = pd.DataFrame({"znox":scale(hprice2["nox"])})
zcrime = pd.DataFrame({"zcrime":scale(hprice2["crime"])})
zrooms = pd.DataFrame({"zrooms":scale(hprice2["rooms"])})
zdist = pd.DataFrame({"zdist":scale(hprice2["dist"])})
zstratio = pd.DataFrame({"zstratio":scale(hprice2["stratio"])})
hprice2 = pd.concat([hprice2, zprice, znox, zcrime, zrooms, zdist, zstratio],axis=1)

housing_standardized = smf.ols("zprice ~ znox + zcrime + zrooms + zdist + zstratio"
                               , data = hprice2).fit()

housing_level.summary()

housing_standardized.summary()

Example 6.2: Effects of Pollution on Housing Prices, Quadratic Interactive Term

Modify the housing model from example 4.5, adding a quadratic term in rooms:

ldist = pd.DataFrame({"ldist":np.log(hprice2["dist"])})
roomssq = pd.DataFrame({"roomssq":hprice2["rooms"]**2})
hprice2 = pd.concat([hprice2, ldist, roomssq],axis=1)

housing_model_4_5 = smf.ols("lprice ~ lnox + ldist + rooms + stratio"
                            , data=hprice2).fit()
housing_model_6_2 = smf.ols("lprice ~ lnox + ldist + rooms + roomssq + stratio"
                            , data=hprice2).fit()

housing_model_4_5.summary()

housing_model_6_2.summary()

Chapter 7: Multiple Regression Analysis with Qualitative Information

Example 7.4: Housing Price Regression, Qualitative Binary variable

This time, use the hrprice1 data.

dataWoo('hprice1', description=True)

name of dataset: hprice1
no of variables: 10
no of observations: 88

+----------+------------------------------+
 | variable | label                        |
+----------+------------------------------+
 | price    | house price, $1000s          |
 | assess   | assessed value, $1000s       |
 | bdrms    | number of bdrms              |
 | lotsize  | size of lot in square feet   |
 | sqrft    | size of house in square feet |
 | colonial | =1 if home is colonial style |
 | lprice   | log(price)                   |
 | lassess  | log(assess                   |
 | llotsize | log(lotsize)                 |
 | lsqrft   | log(sqrft)                   |
+----------+------------------------------+

Collected from the real estate pages of the Boston Globe during 1990.
These are homes that sold in the Boston, MA area.

Estimate the coefficients of the above linear model on the hprice data set.

hprice1 = dataWoo('hprice1')
housing_qualitative = smf.ols("lprice ~ llotsize + lsqrft + bdrms + colonial",
                              data=hprice1).fit()
housing_qualitative.summary()

Chapter 8: Heteroskedasticity

Example 8.9: Determinants of Personal Computer Ownership

Christopher Lemmon, a former MSU undergraduate, collected these data from a survey he took of MSU students in Fall 1994. Load gpa1 and create a new variable combining the fathcoll and mothcoll, into parcoll. This new column indicates if either parent went to college.

dataWoo('gpa1', description=True)

name of dataset: gpa1
no of variables: 29
no of observations: 141

+----------+--------------------------------+
 | variable | label                          |
+----------+--------------------------------+
 | age      | in years                       |
 | soph     | =1 if sophomore                |
 | junior   | =1 if junior                   |
 | senior   | =1 if senior                   |
 | senior5  | =1 if fifth year senior        |
 | male     | =1 if male                     |
 | campus   | =1 if live on campus           |
 | business | =1 if business major           |
 | engineer | =1 if engineering major        |
 | colGPA   | MSU GPA                        |
 | hsGPA    | high school GPA                |
 | ACT      | 'achievement' score            |
 | job19    | =1 if job <= 19 hours          |
 | job20    | =1 if job >= 20 hours          |
 | drive    | =1 if drive to campus          |
 | bike     | =1 if bicycle to campus        |
 | walk     | =1 if walk to campus           |
 | voluntr  | =1 if do volunteer work        |
 | PC       | =1 of pers computer at sch     |
 | greek    | =1 if fraternity or sorority   |
 | car      | =1 if own car                  |
 | siblings | =1 if have siblings            |
 | bgfriend | =1 if boy- or girlfriend       |
 | clubs    | =1 if belong to MSU club       |
 | skipped  | avg lectures missed per week   |
 | alcohol  | avg # days per week drink alc. |
 | gradMI   | =1 if Michigan high school     |
 | fathcoll | =1 if father college grad      |
 | mothcoll | =1 if mother college grad      |
+----------+--------------------------------+

Christopher Lemmon, a former MSU undergraduate, collected these data
from a survey he took of MSU students in Fall 1994.

gpa1 = dataWoo('gpa1')

gpa1['parcoll'] = 0
gpa1.parcoll.loc[gpa1.fathcoll == 1 | gpa1.mothcoll] = 1

GPA_OLS = smf.ols("PC ~ hsGPA + ACT + parcoll", data=gpa1).fit()
GPA_OLS.summary()

Chapter 9: More on Specification and Data Issues

Example 9.8: R&D Intensity and Firm Size

From Businessweek R&D Scoreboard, October 25, 1991. Load the data and estimate the model.

rdchem = dataWoo("rdchem")

all_rdchem = smf.ols("rdintens ~ sales + profmarg", data=rdchem).fit()

all_rdchem.summary()

So, we can estimate the model without that data point to gain a better understanding of how sales and profmarg describe rdintens for most firms. We can use the subset argument of the linear model function to indicate that we only want to estimate the model using data that is less than the highest sales.

rdchem_sub = rdchem[rdchem.sales < max(rdchem.sales)]

smallest_rdchem = smf.ols("rdintens ~ sales + profmarg", data=rdchem_sub).fit()

smallest_rdchem.summary()

Chapter 10: Basic Regression Analysis with Time Series Data

Example 10.2: Effects of Inflation and Deficits on Interest Rates

Data from the Economic Report of the President, 2004, Tables B-64, B-73, and B-79.

dataWoo("intdef", description=True)

name of dataset: intdef
no of variables: 13
no of observations: 56

+----------+----------------------------------+
 | variable | label                            |
+----------+----------------------------------+
 | year     | 1948 to 2003                     |
 | i3       | 3 month T-bill rate              |
 | inf      | CPI inflation rate               |
 | rec      | federal receipts, % GDP          |
 | out      | federal outlays, % GDP           |
 | def      | out - rec                        |
 | i3_1     | i3[_n-1]                         |
 | inf_1    | inf[_n-1]                        |
 | def_1    | def[_n-1]                        |
 | ci3      | i3 - i3_1                        |
 | cinf     | inf - inf_1                      |
 | cdef     | def - def_1                      |
 | y77      | =1 if year >= 1977; change in FY |
+----------+----------------------------------+

Economic Report of the President, 2004, Tables B-64, B-73, and B-79.

intdef = dataWoo("intdef")

intdef['infla'] = intdef['inf']
intdef['defla'] = intdef['def']

tbill_model = smf.ols("i3 ~ infla + defla", data=intdef).fit()
tbill_model.summary()

Example 10.11: Seasonal Effects of Antidumping Filings

C.M. Krupp and P.S. Pollard (1999), Market Responses to Antidumpting Laws: Some Evidence from the U.S. Chemical Industry, Canadian Journal of Economics 29, 199-227. Dr. Krupp kindly provided the data. They are monthly data covering February 1978 through December 1988.

dataWoo("barium", description=True)

name of dataset: barium
no of variables: 31
no of observations: 131

+----------+---------------------------------+
 | variable | label                           |
+----------+---------------------------------+
 | chnimp   | Chinese imports, bar. chl.      |
 | bchlimp  | total imports bar. chl.         |
 | befile6  | =1 for all 6 mos before filing  |
 | affile6  | =1 for all 6 mos after filing   |
 | afdec6   | =1 for all 6 mos after decision |
 | befile12 | =1 all 12 mos before filing     |
 | affile12 | =1 all 12 mos after filing      |
 | afdec12  | =1 all 12 mos after decision    |
 | chempi   | chemical production index       |
 | gas      | gasoline production             |
 | rtwex    | exchange rate index             |
 | spr      | =1 for spring months            |
 | sum      | =1 for summer months            |
 | fall     | =1 for fall months              |
 | lchnimp  | log(chnimp)                     |
 | lgas     | log(gas)                        |
 | lrtwex   | log(rtwex)                      |
 | lchempi  | log(chempi)                     |
 | t        | time trend                      |
 | feb      | =1 if month is feb              |
 | mar      | =1 if month is march            |
 | apr      |                                 |
 | may      |                                 |
 | jun      |                                 |
 | jul      |                                 |
 | aug      |                                 |
 | sep      |                                 |
 | oct      |                                 |
 | nov      |                                 |
 | dec      |                                 |
 | percchn  | % imports from china            |
+----------+---------------------------------+

C.M. Krupp and P.S. Pollard (1999), "Market Responses to Antidumpting
Laws: Some Evidence from the U.S. Chemical Industry," Canadian Journal
of Economics 29, 199-227. Dr. Krupp kindly provided the data. They are
monthly data covering February 1978 through December 1988.

barium = dataWoo("barium")

barium_imports = smf.ols(
    "lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6",
    data=barium).fit()
barium_imports.summary()

Estimate a new model, barium_seasonal which accounts for seasonality by adding dummy variables contained in the data.

barium_seasonal = smf.ols(
    "lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6 + \
    feb + mar + apr + may + jun + jul + aug + sep + oct + nov + dec",
    data=barium).fit()
barium_seasonal.summary()

Now, compute the anova between the two models.

from statsmodels.stats.api import anova_lm

anova_lm(barium_imports, barium_seasonal)

in conclusion

--I would be grateful if you could point out or comment on Qiita or github. -I also do Twitter.