Theoretical foundations of industrial safety in modern conditions

The service sector in the industrialized countries throughout the world has taken a leading position in the socio-economic life. Thereby, ensuring the improvement of the provided service quality and the expansion of service range are becoming a nationwide business. Service is the economy development as well as the public employment and the basis of stability in society. Moreover, service should be under control; and one of the most efficient ways to provide it with supervision is an expert review. In all countries with a market economy, such expert review (the notion of “inspection” is used abroad) has long taken its place and its influence on the regulation of production activities is quite favorable [1].

Introduction
The service sector gives from 60 to 75% of the national income and employs more than 70% of the working-age population. In Russia, this area provides about 1/2 of the national income. When joining the WTO, domestic enterprises will be obliged to withstand strong competition from Western partners as well as the expansion of foreign goods in the markets.
Due to the rapid development of market relations on the territory of the former USSR, the liberalization of foreign trade and the weakening of the administrative regulatory state role, the traditional market regulation and control methods for product reproduction are becoming more important. One of the most important mechanisms for such regulation and control is an objective and competent expert review. The world-famous companies such as "Lloid", "Veritas", TUV and others, as well as the associations of inspection companies, for example "IFIA", operate in this area.

Materials and methods
As for the Russian industry, the use of the expert review instrument is carried out within the framework of the industrial safety management system. The industrial safety expert review, in accordance with the law "On industrial safety of hazardous production facilities", forms the basis of industrial safety, being an integral safety management instrument. For the first time in Russia, the concept of "industrial safety expert review" was introduced by the Federal Law of July 21, 1997 No. 116-FL [2], expanding the previously existing concept of the expert review of industrial plants (facilities), equipment and works that were under the jurisdiction of the activities licensed by the Russian Federation Oversight Committee for Mining and Industry. The main expert review tasks include determining the safety degree of the decisions in the area of design and all other documentation regarding technical devices, buildings and structures at a hazardous production facility; and also the industrial safety expert review declarations and other documents related to the operation of a hazardous production facility.
The industrial safety expert review is carried out in accordance with the rules established by the Federal Service for Technological, Environmental and Nuclear Supervision [3] by organizations that have licenses to carry out industrial safety expert review at the expense of the organization (customer), assuming or operating a dangerous production facility. The industrial safety expert review declarations is carried out in accordance with the rules [4,5].
The most important condition for the examination is its reliability, independence and enhanced quality of expert work. The fulfillment of these conditions is achieved by the proper organization of expert work, the selection of qualified personnel, the provision of expert organization of normativetechnical and legal documentation, technical equipment, awareness of economic and legal responsibility. The regulation of activities in the industrial safety expert review by the state is carried out through the scope of licensing (issuance of documents authorizing this activity type) and subsequent supervision (control) over compliance with the requirements and conditions by the license holder. Such regulation does not completely solve the problems arising in the quality of the industrial safety expert review.
The law "On Technical Regulation" provides for the concept of accreditation that is "an official recognition by an accreditation body of the competence of a natural or legal person to perform work in a specific area of conformity assessment. Conformity assessment is a direct or indirect determination of compliance with the requirements for the object" [6].
The Supervisory Board is formed from representatives of the Federal Service for Environmental, Technological and Nuclear Supervision (FSETNS), its territorial bodies and its subordinate organizations. It monitors the activities of the commission of ISER. The composition of the Supervisory Board is approved by FSETNS. The Advisory Board consists of representatives of organizations interested in the work of the commission of ISER and has a consultative view. The composition of the Advisory Board is approved by the Supervisory Board. Sectoral commissions are created by the Supervisory Board and solve specific, professional tasks in areas corresponding to their competence. The coordinating body ensures the detailing of the entire Accreditation System, in particular, it analyzes the performance of expert organizations, the regulatory and methodological basis of the commission of ISER.

Results and discussion
The industrial safety expert review is carried out at the expense of the customer on the basis of his application or other documents in accordance with the conditions agreed by the expert organization and the customer. The expert organization starts to conduct an inspection only after receiving a set of necessary materials and documents in full in accordance with the requirements of the current regulatory technical documents. The result of the industrial safety expert review is the conclusion. The decision to issue a positive or negative expert opinion is made on the basis of the review and analysis of the documents obtained during the expert review, checking the condition of the object or conducting the necessary tests.
The organizations involved in the provision of industrial safety expert review services are expert organizations that are business services enterprises. The creation of such service enterprises is an integral trend of a market economy. In today's economic environment, consumer requirements, supported by relevant legislation on consumer protection, are becoming increasingly weighty and stringent, which leads to significant changes in the structure of product quality indicators, and, consequently, to an increase in the socio-economic significance of the problem of improving the quality of technical services. Experts point to different attitudes of consumers to assess the quality of products and to assess the quality of services. The consumer considers the quality of the service not only as a result of the receipt of a product or service, but also as the process of providing this service  3 and the nature of the interaction that takes place. Thereby, quality management services is a connecting element of each process of activity and an integral component of the overall management of the technical service enterprises. Wherever possible, the process, product and service should be certified.
The degree of implementation and the nature of the production goals set by the enterprise help to make economic and scientific-technical forecasts. Forecasting allows you to identify ways and directions of development of the enterprise, to determine the timing of their implementation. Initially, a technical forecast is carried out, the specific needs of the market are identified, in this case, in the production of safe products and meeting certain industrial safety requirements. The main operational parameters of technical devices are determined in accordance with the requirements of the industrial safety expertise conducted at the design stage. Economic forecasting is carried out in order to identify the needs of the market for a new product; the cost of manufacturing the product and the price on the market (taking into account the costs of meeting the requirements of technical regulations and conducting expertise). Thus, the technical and economic forecast in this case is made on the basis of the specific needs of the market for technical devices used in the HIF, describes possible methods (technologies, organization) of production and the development of the object of forecasting.

Conclusion
Since the marketing as a creative management function intends to conduct scientific research aimed at meeting the needs and requirements of customers, the results obtained at the forecasting stage using industrial safety expert review and the study of the market requirements (state, legislation) for the manufactured goods (equipment, technical devices ), are an essential element of marketing research. At the stage of operational marketing, manufacturers determine the ability of the manufactured equipment to undergo industrial safety expertise. The required characteristics of safe operation are determined at this stage of the product life cycle.
Certain safety requirements are also imposed on the material and technical supply of the production of machinery in accordance with the law on technical regulation. Such requirements include the availability of a certificate of conformity from suppliers of raw materials, i.e. Confirmation of compliance of objects with the requirements of technical regulations ensuring industrial safety.
At the design stage of the production process (preparation and development of the production process), among other things, the tasks of ensuring the established profile and characteristics of the manufactured product and compliance with technical regulations are solved. Thus, the safety expert review, carried out at the design stage, and high-quality, certified raw materials obtained at the material and technical supply stage shall make it possible to develop an industrial process that is safe for society and the environment.
At the next phase of the product life cycle, at the stage of production of technical devices, the Law "On Technical Regulation" establishes requirements for the processes of production of technical devices used in HIFs. It comes to the established special technical regulations (art. 8) for those particular types of products and their processes of production, operation, storage, transportation, sale and disposal, the risk of harm being higher than the risk of harm taken into account by the general technical regulation. This type of product includes technical devices used in HIFs (Table 1).
It has been revealed that at all stages of the life cycle of technical devices used at HIFs, within the framework of the legislation of the Russian Federation, ISER is necessary. The expert review is one of the instruments to ensure industrial safety, along with standardization, certification and declaration of industrial safety. The industrial safety expert review allows the state to control the level of risk for employees of the enterprise (reduce it) and the degree of environmental risk in the production of technical devices.
Certainly, for the manufacturer of the equipment for obtaining certificates of conformity, fulfilling the requirements of technical regulations and forecasting the requirements imposed within the ISER, there is an increase in production costs (increased costs for ensuring industrial safety, advertising, creating a brand of such a high-level manufacturing company, a distribution network and etc.), but also:  increase the efficiency of production activities. Due to the rising costs of ensuring industrial safety, the quality of production of technical devices improves, this leads to an increase in the price of products; for the consumer (the company operating the HIF and the technical devices used on it), working with a manufacturer of this level allows you to minimize the risk and possible damage from the consequences of accidents;  advertising and brand awareness for the manufacturer contribute to an increase in market share, sales volume and level of competitiveness in the market for goods. Accounting of possible risks and accidents in the declaration of industrial safety subject to industrial safety expert review The Law on Industrial Safety also establishes a number of requirements for conducting industrial safety expertise and defines the procedure for its implementation, i.e., in fact, the requirements for the provision of services are established. A number of specialists consider the industrial safety expert review of technical devices as an assessment of product conformity [7,8]. In this case, the timing, procedure for conducting the safety expert review itself are also contained in the technical regulations.


The Khaldun-Laffer Curve Revisited: A Personal Income Tax-Based Analysis for Turkey

The objective of this paper is to revisit as well as empirically examine an old but still discussed postulate, the Khaldun-Laffer curve, on the basis of personal income tax by making use annual time-series data for Turkey for the period 1970-2015. The findings of the paper confirm the validity of the Khaldun-Laffer curve hypothesis. In addition, we infer that the optimal tax rate that maximizes the tax revenue generated from personal income taxation in Turkey is 15.03 percent. This rate is well-below than the current rate which we estimate as 15.37 percent, implying that Turkey’s current tax rate for personal income tax takes place in the prohibitive range of the Khaldun-Laffer curve. These findings suggest that the current tax rate should be lowered and to its optimal level to collect more tax revenue. Getting down the current rate to its revenue-maximizing rate not only would it enable the Turkish authorities to collect more revenues with a relatively lower rate, but also would allow them to minimize the substitution effects of personal income tax while maximizing the income revenues from it.

Introduction
One of the most important and well-known arguments of supply-siders is the notion of the Khaldûn-Laffer curve (henceforth the K-L curve) 1 ; based on persistent high tax rates in the US, supply-siders have claimed that lower tax rates would mean higher tax revenue and that consequently higher rates would result in lower tax revenue. Hence, reducing current tax rates would be a good tax policy because it will not only maximize tax revenue but also motivate economic activity. This argument of supply-siders made the K-L curve popular in the early 1980s. The K-L curve is the graphical representation of the aforementioned postulate. It depicts a classic bellshaped relationship 2 between tax rates and tax revenue: as the tax rate increases the tax revenue increases till it reaches an optimal point of maximum tax revenue; any further increase in the tax rate decreases the revenue due to the disincentive effects created by higher taxes. Laffer (1981Laffer ( , 2004 argues that higher tax rates will discourage work and production by removing incentives, triggering sluggish growth, and thereby lowering tax revenue collected by the government. As with many other developing countries, maximizing revenue is important for Turkey in order to realize its fiscal targets. Considering the undesired consequences of the alternative ways to finance budget deficits (especially, borrowing and printing money), the significance of having an ideal tax system in which tax rates generate maximum tax revenue for the government without inducing a dampening effect on the economy has been well understood.
Turning to the case of Turkey, this paper proposes to revisit the K-L curve. In this context, the paper first endeavors to examine the validity the K-L curve empirically and then attempts to find the optimal tax rate at which the maximum amount of tax revenue can be collected. The paper also aims to: (i) investigate whether present tax rates are lower or higher than the revenue-maximizing tax rates; (ii) determine the income and substitution effects of tax rates and thereby design ideal tax rates for fiscal policy purposes; (iii) avoid the undesired consequences of taxation incurred from high tax rates (i.e. tax avoidance, tax evasion, preferring leisure to work, changing the scope and/or structure of economic activity, etc.) which result not only in a reduction of tax revenue but also negatively affect economic activity and thus growth, and (iv) contribute to the limited empirical literature related to the K-L curve for developing countries like Turkey.
Since the curve theoretically establishes a non-linear link between tax rates and tax revenue, we have chosen the personal income tax (hereafter, the PIT) as the case for testing the validity of it. The reasons for doing so, first and foremost, is that the PIT is a tax that is levied on individuals' taxable income. Second, this tax accounts for the largest share of the central government's tax revenue after VAT and the special consumption tax. Thirdly, the PIT is a tax instrument that is well-suited to examining the K-L curve.
The PIT is one of four important taxes that generate the highest tax revenue in the Turkish tax system. The other taxes are the corporate income tax, the value-added tax, and the special consumption tax. These four taxes account for an average of more than 80% of total central government tax revenue. Although two indirect taxes (the value-added tax and special consumption tax) dominate the tax system and constitute nearly 65-70% of the central government's tax revenue on average, the PIT also has a significant place in the tax system, accounting for almost 20% of it. Moreover, it is the sole progressive tax (if wealth tax, which constitutes a minor portion of total tax revenue, is ignored). At present, the PIT contains four brackets with rates of 15, 20, 27 and 35%, respectively. As a rule, these rates are not frequently changed. The latest alteration was made in 2006 when the number of tax brackets was reduced from 6 to 4. However, at the end of the year, before starting the new fiscal year, tax brackets are re-adjusted for subsequent fiscal years in line with developments in inflation of the current year.
The remainder of the paper is organized as follows. Section 2 outlines the theoretical background and related literature. Section 3 describes the data and econometric methodology, while Section 4 presents the estimation strategy and results of the paper alongside their interpretation; the final section concludes.

Theoretical background
While Laffer was discussing the US President's proposal for tax increases in a dinner meeting in 1974, he charted a classic bell-shaped curve on a cocktail napkin, illustrating the trade-off between tax rates and tax revenue (see Laffer, 1981Laffer, , 2004Wanniski, 1978Wanniski, , 2005. Soon after, Wanniski (1978) named this trade-off as 'the Laffer curve'; however, Laffer (1981Laffer ( , 2004 notes that the idea behind the curve belongs to Ibn Khaldûn, a 14th-century Muslim philosopher. Giving a special reference to the opus of Khaldûn, the Muqaddimah, Laffer draws attention to the background of the curve with the words of Khaldûn himself: 'it should be known that at the beginning of the dynasty, taxation yields a large revenue from small assessments. At the end of the dynasty, taxation yields a small revenue from large assessments' (Khaldûn, 1980 [1377], p. 80).
However, contrary to the arguments above, some authors, such as Fullerton (1982), van Ravestein and Vijlbrief (1988), and Hsing (1996), attribute the roots of the curve to Adam Smith, an 18th-century Scottish economist, referencing his words: 'high taxes, sometimes by diminishing the consumption of the taxed commodities, and sometimes by encouraging smuggling, frequently afford a smaller revenue to government than what might be drawn from more moderate taxes' (Smith, 2008(Smith, [1776, p. 835).
Despite these counter-arguments, we prefer to label the curve as the K-L curve to give credit to both Khaldûn and Laffer.
The curve exhibits a one-way causality from tax rates to tax revenue, illustrating a hump-shaped relationship between tax rates and tax revenue. The upward-sloping side of the K-L curve displays the normal range, whereas the downward-sloping side shows the prohibitive range (Wanniski, 1978;Laffer, 1981Laffer, , 2004. This implies that a given tax revenue can be collected with two different tax rates, one falls into the normal range of the curve, and the other one falls into the prohibitive range; the rational fiscal policy option would be choosing the lower rate, which falls into the normal range because the income/incentive effect dominates the substitution/disincentive effect of taxation.
The K-L curve is still being discussed, although more than three decades have passed since its postulation by Arthur Laffer. The focal point of discussion is whether the K-L curve is true or false, or only an approximation. The discussions are centered specifically on whether there is such a relation between tax rates and tax revenue, as emphatically argued by Laffer (1981Laffer ( , 2004. If so, what is the revenue-maximizing tax rate, i.e. the optimal tax rate? The K-L curve received a great attention especially in the early 1980s with the arguments of the prominent supply-siders such as Paul Evans, Michael Boskin and Martin Feldstein, asserting that 'lower tax rates would mean higher revenue because existing rates were too high to maximize tax revenue, -that is, tax rates were so high that fewer taxed goods were being produced and the overall effect was lower tax revenue' (Becsi, 2000, p. 53). Even today, the curve still receives remarkable attention in tax-related discussions, especially in the United Stated of America, where a conservative political movement within the Republican Party that is called 'Tea Party' keeps the topic always on the agenda and makes the reduction of taxes one of its primary objectives (Tanzi, 2014).
Despite these facts, the curve encounters severe criticism in the literature, such as: (i) the K-L curve, in the words of Tanzi (2014), is a wonderful propaganda device for conservatives; (ii) the idea that the Laffer curve has disincentive effects on high tax rates is not precise, and (iii) there are suspicions about the existence of a relationship between tax rate and tax revenue. More specifically, the curve that establishes a classic bell-shaped relationship between tax rates and tax revenue obtained by the government may not reflect the reality. Depending on a number of factors, such as taxpayers' reaction to the taxes, the openness degree of the economy, efforts of the tax administration, structure of labor market, and etc., the curve may take different shapes, rather than proposed by Arthur Laffer 3 .
In the literature, we see at first glance that there have been many studies on the K-L curve, a large majority of which have focused on the theoretical aspects of the curve, rather than empirically exploring it. There are also numerous counter-studies such as Henderson (1981), Sanyal, Gang and Goswani (2000), Busato and Chiarini (2013), which critically assess the curve in terms of different aspects, ranging from its assumptions to its validity.
This paper contributes to the literature not only by focusing on the PIT of Turkey in the context of the K-L curve in the long run, but also by including control variables other than tax rate and its square as determinants of tax revenue in order to prevent omitted variable bias that may produce misleading results.

Related literature
An early study by Stuart (1981) calibrates the Laffer curve for Sweden for the 1970s and then concludes that the Swedish tax rate on labor income is higher than the revenue-maximizing tax rate. He estimates the revenue-maximizing marginal tax rate between 69-73%, while the actual tax rate is around 80%. These findings confirm that the Swedish tax system was in the prohibitive range of the K-L curve in the 1970s. On the other hand, Feige and McGee (1983) find that the marginal tax rate was about 0.83, whereas the revenue-maximizing rate was 0.58 for Sweden for 1979. Both studies show that the optimal tax rate is lower than the actual one for Sweden. This implies the importance of cutting tax rates for this country for obtaining a revenue-maximizing rate as well as reducing the substitution effects of taxation arising from high tax rates. Van Ravestein and Vijlbrief (1988) employ a similar model to Stuart (1981) and then estimate the actual tax rate and the peak level of the Laffer curve for the Netherlands; they found that all the actual tax rates are lower than the revenue-maximizing rate, reflecting the availability of room for hiking the marginal tax rates for this country.
Another study by Hsing (1996) investigates the Laffer curve for the US by employing a single-factor model, considering tax revenue as a function of only tax rate and its square; the author finds evidence that the revenue-maximizing tax rate changes for the PIT between 32.67% and 35.21% for this country for the period 1959-1991. On the other hand, the average tax rate is 19.58% based on tax liability and 20.18% based on tax payment in 1991. Based on these results, Hsing (1996) argues that there is a room for increasing the average tax rate to collect more tax revenue. Similarly, Karas (2012) examines the revenue-maximizing PIT rate for the Czech Republic. Karas (2012) confirms the existence of an inverted U-shape relationship between the tax rates and tax revenue; the revenue-maximizing PIT rate for the Czech Republic is 33.13%. Walewski (1999) explores the existence of the K-L curve by considering the case of three European transition economies (the Czech Republic, Poland, and Hungary) and observes that the curve for the selected countries can be classic bell-shaped. Based on the findings, both the Czech Republic and Poland are still on the upward sloping side of the curve, whereas Hungary swings around the Laffer hill: that refers to the revenue-maximizing point of the curve.
The existing literature contains also some other prominent studies that focus on country groups. For instance, Dalamagas (1998) examines the Laffer curve for 13 industrialized OECD countries and finds that the curve has the familiar inverted U-shape; his findings also document that a permanent reduction in tax rates has a potential of generating more tax revenue in the long run for the highly taxed countries alone. However, he argues that a permanent reduction in average tax rates may result in increases in the long-run of government budget deficits in economies with crowding-out characteristics, thus providing a strong refutation of the dynamic Laffer curve proposition.
Two follow-up studies by Uhlig (2011, 2013) investigate quantitatively the shape of the Laffer curves for labor taxation and capital income taxation for a group of countries covering the US, EU-14 and some other non-EU countries. The aforementioned authors document that the Laffer curve exists for labor taxes as well as capital taxes. According to Uhlig (2011, 2013), both the US and EU-14 are located on the left side of the curve, but the EU-14 countries are much closer to the slippery slopes than the US. More specifically, they suggest that increasing tax revenue by 30% in the US requires raising labor taxes but only 6% by raising capital income taxes. However, the same requirement for the EU-14 is 8% and 1%, respectively. Another interesting finding of their study is that in the US, 32% of a labor tax cut and 51% of a capital tax cut are self-financing in the steady state, whereas in the EU-14 countries, 54% of a labor tax cut and 79% of a capital tax cut are self-financing.
Likewise, Nutahara (2015) explores Laffer curves for labor, capital and consumption taxes for Japan and his findings confirm that the Laffer curves for labor and capital taxes have single peaks, whereas consumption tax revenue increases monotonically with respect to the tax rate. Furthermore, according to the findings of the author, while the labor tax rate is smaller than that at the peak of the curve, the capital tax rate is either very close to, or bigger than, that at the peak of it.
In sum, a large number of empirical studies on the K-L curve concentrate on either an industrialized country or on country groups, while studies that focus on the K-L curve for other countries, i.e. developing, emerging, and transitional, are highly limited, including Turkey. To the best of our knowledge, there are few empirical works that examine the Laffer curve for Turkey. However, in contrast to the present paper, all the available studies without exception, do not take the issue on a specific tax as we did and employ mostly univariate models.
As stated above, with few exceptions, virtually in all studies different models have been used and this makes it difficult to make comparison among countries, unless the same models are used for different country groups. On the other hand, most of the assumptions of the models employed in the empirical studies are highly far from the reality. In other words, they do not match with the actual economy and taxation system.
Contrary to many other studies, in this paper, we specifically examine the K-L curve in the context of the PIT. Secondly, as opposed to a large strand of studies that only consider tax rate and its square as independent variables for tax revenue we employ a long-run estimation analysis with multiple variables (the share of trade, agriculture, external debt in GDP, inflation rate and GDP per capita) to avoid misspecification problem. One more merit of our paper is that we work with relatively longer time span with 45 annual observations.

Data
The data used in this paper spans the 1970-2015 period 4 ; the availability of the data is the main reason for choosing the study period. The definitions and the sources of the variables used in our model are presented in Table 1. Increasing the precision of the econometric methodology highly depends on the number of observations, so we try to keep the time horizon as large as possible to increase the reliability of our econometric model. Large sample size leads to increased precision in estimates of population. So, 'more is better' phrase suits well. In addition, forecasting by using time series analysis requires higher sample size. However, in this paper we do not focus on forecasting. Rather, concentrating on finding out an optimal tax rate. For many developing countries, lack of data is observed excessively due to several reasons, such as making reforms more frequently, changing the calculation methods in data collection or collecting data for different time lines. As a result, analysing for a developing country case is challenging. In our analysis, the number of observations and degrees of freedom are larger than 30. We believe that statistically all these are fair enough for interpreting the coefficients.

Econometric methodology
In this section, besides the econometric methodology for our model, we also briefly provide the theoretical framework for the variables we chose. As noted earlier, there are several theoretical formulations for the K-L curve in the literature, each of which is based on different assumptions. Existing empirical studies on the curve typically employ a single factor model 5 . However, omitting a relevant variable causes specification error, the so-called 'omitted variable bias'. Underspecifying the model by excluding a relevant variable violates one of the Gauss-Markov assumptions: zero conditional means of the error term (Woolridge, 2012). This violation harms the unbiasedness of estimated coefficients. Hence, if the expected values of estimated coefficients are not equal to population parameters, it will be meaningless to interpret them; in order to avoid this misspecification problem, some control variables are included in our analysis.
In the present paper, we depict the K-L curve by controlling other determinants of PIT revenue different from its tax rate. These control variables are as follows: (i) GDP per capita; (ii) share of agriculture; (iii) share of trade; (iv) share of external debt; and (v) the inflation rate. With the exception of GDP per capita and the inflation rate, all other control variables are expressed as a fraction of GDP.
Since GDP per capita is simply the overall production per person, it can be used as a proxy for the whole economy (Gupta, 2007). Then, we include inflation rate in our analysis to capture the price effects in the economy; for this purpose, we define the inflation rate as the annual percentage change in the consumer price index over the previous year. Yet another control variable is the share of agriculture in GDP, which represents the informal sector that cannot be taxed properly from personal income, owing to lack of well-functioning bookkeeping system and market economy, dominance of small-sized enterprises, and existence of subsistence farming make difficult to tax the agriculture sector for a large number of developing countries. Furthermore, in parallel with the extant literature, the share of trade in GDP is used as a proxy of openness of the economy to the international trade. The last control variable, external debt share, indicates the government's revenue generating performance (Gupta, 2007).
The challenging point in this paper is to define the tax rate since the Turkish PIT has a progressive tax structure; depending on the level of the taxpayers' current income, the tax rates are set at 15, 20, 27 and 35%. Since it is practically almost impossible to access the number of taxpayers who fall in which tax bracket during the fiscal year, we choose using the PIT revenue for each employed -to -GDP per capita ratio which can be considered as the individual's income-tax burden. This kind of calculation methodology is frequently applied in the literature, especially in the absence of a specific flat tax rate as in the case of Turkey. The calculated value can be called the marginal tax rate. We use PIT revenue as the dependent variable and to eliminate the price effect, we converted the nominal level of tax revenue into the real one.
We aim to find the threshold level for the tax rate when the behavior of tax revenue changes by controlling the following variables. The appropriate econometric model for this purpose can be designed as follows: where T and tax denote PIT revenue and rate respectively, whereas X denotes control variables which are listed before so that b 0 , b 1 , b 2 are scalars and b 3 is a parameter vector. Proving the possible existence of the curve for Turkey, the expected sign of the coefficient on PIT rate (b 1 ) should be positive while the coefficient on PIT rate square (b 2 ) should be negative. Then, the optimal tax rate will be:

Estimation strategy and results
The problem of spurious regression may arise when time series data is applied to the non-stationary form; in such a case, one possible solution is to make the series stationary by taking first-differences. However, it is important to keep in mind that the differencing of the series would hinder long-run analysis (Davidson et al., 1978). To avoid such a problem, there are several methods. For instance, Davidson et al. (1978) propose an error correction mechanism by using first differences of the short-run and undifferenced values for the long-run dynamics. However, Engle and Granger (1987) prove the necessity of cointegration relations in order to be able to implement this method. The Johansen Cointegration test is preferred to realize the long-run relationship between variables. To do so, the first step is to test the variables in terms of stationarity.

Unit root analysis
The null hypothesis of the Dickey-Fuller (DF) unit root test states that the data needs differencing to be stationary, whereas the alternative hypothesis argues the opposite; in other words, data may not require differencing to be stationary. According to the results of the ADF unit root test, we cannot reject the null hypothesis that implies all variables are non-stationary. To see the integration order, ADF unit root test is applied to the first differences; Table 2 reports the results of this test indicating that all the differenced variables are stationary. Since all the variables are I(1), it is suitable to further apply cointegration analysis.

Cointegration tests and results
The cointegration test is based on the idea that if a linear combination of non-stationary (I(1)) variables leads to stationarity (I(0)), then the variables are cointegrated. Cointegration relation between these variables implies the existence of a long-run equilibrium relationship between them. Johansen (1988), and Johansen and Juselis (1990) construct an unrestricted VAR model by deriving maximum likelihood estimators of cointegrating vectors for an autoregressive process in order to find the number of cointegrating relations by using error correction model 6 .
One of the most crucial assumptions of the Johansen Cointegration test, that all variables are in the same integration order, is satisfied (see Table 2). The Johansen Trace test results confirm that there is a strong long-run relationship between the variables in our model (see Table 3) 7 . Table 3 also tabulates three different Johansen Trace test results. In the first one, all independent variables are included, whereas in the second one external debt share is excluded. In the last one, additional to external debt share, the inflation rate is omitted 8 . Note: ** There is one cointegrating relationship at the 95% confidence level.

Source: Authors' calculations
The findings of our long-run estimates presented in Table 4 are threefold. According to the first regression, the share of agriculture, trade and external debt are insignificant and all other variables are significant at the 95% confidence level (Table 4). In the second regression, the external debt share is omitted since it is statistically insignificant. In the last regression, we exclude the inflation rate with the same concern as in the external debt share. Note: Robust standard errors are in parentheses. a Optimal tax rate level is the revenue-maximizing rate of the K-L curve. b The observation period is 1970-2015. However, one lag is imposed in the vector error correction model by AIC and SBIC criteria. Also, standard deviations are in parentheses. *** p < 0.01, ** p < 0.05.

Source: Authors' calculations
Before turning to discuss the optimal PIT rate, it would be useful to interpret the signs of the control variables included in the models. Wagner and Weber (1977) declare that as an economy progresses, it is also observed that the public sector tends to expand. In the literature, this principle is called Wagner's Law; a positive and significant relationship between PIT revenue and GDP per capita supports the Wagner's law since GDP per capita is used to represent the condition of the overall economy.
On the other hand, the sectorial decomposition of the output is significant for tax revenue because it renders possible to collect further taxes from certain industries (Stotsky and WoldeMariam, 1997;Gupta, 2007). Deficiencies in bookkeeping makes the agriculture sector hard to collect tax properly, especially in developing countries. The significant and negative sign in the last two regressions for agriculture share of GDP points out the adverse relationship between the agricultural sector and PIT, which is consistent with the findings of Stotsky and WoldeMariam (1997) and Gupta (2007). In essence, as the agriculture share of Turkey increases, the Turkish government's revenue performance for personal income taxation declines.
The relationship between the inflation rate and PIT revenue is straightforward since high inflation rates imply rises in the wages in nominal terms. In the view of the progressivity of PIT, inflation drives incomes up and into higher tax brackets, resulting in bracket creep depending on the inflation's speed and the structure of tax progressiveness. Hence, the bracket creep allows the government to generate artificially more PIT revenue by boosting the real tax burden of personal income taxpayers, a phenomenon called 'taxflation' in the literature. On the other hand, it is also very likely that inflation may reduce the real value of tax revenues by eroding tax base, the so-called 'Tanzi effect' or 'the Olivera-Tanzi effect' that emerges in the existence of the following three conditions: (i) high inflation; (ii) lags in tax collection in un-indexed systems, and (iii) inelastic tax structure 9 . Inflation is also likely to erode the tax base by shrinking the formal sector at the expense of the informal sector. Briefly, in theory, inflation may positively or adversely affect PIT revenue depending on the case we explained above. However, our analysis demonstrates that the adverse behaviors of taxflation and the Olivera-Tanzi effect on PIT revenue make the effect of inflation on the revenue statistically insignificant for Turkey. Gupta (2007) argues that when the government needs further revenue to finance its expenditures and/or to pay its due debt, the government will need more tax revenue, and external debt share shows this relationship. In our analysis, we find that the relationship between the external debt share and PIT revenue is too weak, i.e. insignificant. So, it would not be wrong to say that the Turkish government does not use a significant amount of PIT revenue to pay its external debt.
Lastly, as known, international trade is likely to affect tax revenue through various channels. The significant and positive sign of trade share in GDP shows consistency with these channels for Turkey. Firstly, the positive effect of international trade on economic development is a commonly regarded idea since the time of Adam Smith; accordingly, as a country's income rises, revenue obtained from PIT increases, too. Secondly, compared to domestic activities, taxing international trade is relatively easier since international activities are taking place at specific locations of the country.
Our long-run estimates from the last regression 10 show that the revenue-maximizing PIT rate (the optimal PIT rate) is 15.03%. The actual tax rate, on the other hand, for the year 2015 is estimated at 15.37%. Based on these findings, we can argue that Turkey is on the right side of the K-L curve referring to the prohibitive range of the curve (see Figure 1).
This finding suggests that if the Turkish fiscal authorities should aim to increase the PIT rate further with the purpose of generating more revenue from this tax, this tax policy action will end up with a significant reduction in PIT revenue rather than a rise; in order to avoid this situation, the current PIT rate should instead be lowered to the rate that corresponds to the Laffer hill, which is 15.03%.

Source: Authors' computation
10 From now on, we consider the results of the last regression since it performs along the lines of our suggested theoretical framework. 11 Based on Table 3, we construct Figure 1 depending on the coefficients coming from long-run estimates, under the assumption of zero control variables since the aim is only to see the classic bell-shaped relationship between tax rate and tax revenue.

Conclusion
This paper revisited and then empirically tested the potential existence of the K-L curve for personal income tax (PIT) on the long run by utilizing annual Turkish data over the 1970-2015 period. Our estimation results appear to confirm the validity of the K-L curve in the context of the PIT for Turkey; we found the optimal PIT rate as 15.03% and the present PIT rate as 15.37% for Turkey. Thus, we can argue that the current PIT rate in Turkey is somewhat over its optimal, putting Turkey's current PIT rate into the prohibitive range of the curve. This implies that a reduction in the PIT rate increases PIT revenue that will be collected by the Turkish authorities; in other words, if they reduce the current PIT rate down to its optimal, ceteris paribus, PIT revenue will rise to the corresponding maximum level. If done so, the authorities would have an opportunity of generating further tax revenues with a lower tax rate as higher tax rates create a substitution effect that would discourage working as well as forcing taxpayers to seek out the ways of reducing the tax burden. Our findings presented above are broadly consistent with the findings of similar but scant, empirical studies in this field, such as Hsing (1996) and Karas (2012), which solely focus on the PIT.
Realizing what drives tax revenue is crucially important for not only designing but also implementing a proper tax policy for all countries. A good starting point is the K-L curve that relates tax rates with tax revenue; for the PIT, the curve represents a concave relationship between the revenue obtained from this tax and its rates. From a tax policy perspective, on the basis of our empirical findings, the closeness of the current PIT rate (15.37%) to the lowest level of tax bracket that is 15% as highlighted earlier suggests that the main source of the PIT revenue comes from individuals at the lower income level. As such, individuals with lower income whose taxable incomes fall in the first two brackets of the PIT largely bear the tax burden of the PIT, which may explain the distorted redistributive role of the PIT in Turkey. Besides, the optimal PIT rate also implies that in order to maximize tax revenue, the government should not exceed the lowest level of the income bracket.
Despite the positive explanations in favor of the K-L curve throughout the present paper, as frequently discussed elsewhere in the literature, it is also obvious that the curve has several drawbacks and deficiencies on both theoretical and empirical grounds that should be considered for future studies. First, in Laffer's view, the shape of the curve is solely determined by tax rates. However, the curve may take different shapes and positions depending on various factors, such as the state of business cycle, openness to trade, reaction of taxpayers to the changes in tax rates, how the government spends its tax revenue (investment, consumption, or transfers?), tax consciousness and ethics, and the labor market structure and arrangements, and the size of informal economy. The factors above and their similarities can also make the curve unstable, changing the slope of the curve subject to the time horizon. Second, it is uncertain whether cuts in tax rates produce greater tax revenue while boosting economic activity by incentivizing saving and investment as suggested by Laffer. Third, Laffer does not mention whether tax rate reductions are permanent or transitory. We know that the effects of each measure on the economy are different. Fourth, tax rate cuts would not always encourage the taxpayer to work more as proposed by Laffer. Last, there may be more than one peak of the K-L curve that maximizes tax revenue.
Having considered the aforementioned drawbacks and deficiencies, the Laffer curve is still a critically important concept that draws our attention to the harmful effects of high taxes not just for revenue enhancing, but also for stimulating economic activity by introducing an optimal tax rate. No doubt, further empirical research is essential to identify the drawbacks and deficiencies of the curve for having better tax policies in particular and for designing growth-friendly economic policies in general.


Causal inference of hazard ratio based on propensity score matching

Propensity score matching is commonly used to draw causal inference from observational survival data. However, there is no gold standard approach to analyze survival data after propensity score matching, and variance estimation after matching is open to debate. We derive the statistical properties of the propensity score matching estimator of the marginal causal hazard ratio based on matching with replacement and a fixed number of matches. We also propose a double-resampling technique for variance estimation that takes into account the uncertainty due to propensity score estimation prior to matching.

Introduction
For observational survival data, Cox proportional hazards modeling is the most common approach to infer the association between treatment and survival outcome (Cox;1972). Following Rubin's causal model, Hernán et al. (2000) and Robins et al. (2000) proposed the marginal structural Cox proportional hazards model and inverse probability of treatment weighting to estimate the causal effect of treatment on survival outcome.
Under the classical assumptions of no unmeasured confounding and positivity, the causal effect of treatment is identifiable from the observed data. The seminal paper Rosenbaum and Rubin (1983a) has demonstrated the central role of propensity score as the balancing score, in the sense that balancing the propensity score can balance the observed covariates on average. Researchers have proposed regression, propensity score-based weighting, stratification, or matching estimators to adjust for confounding. In the context of survival analyses, Chen and Tsiatis (2001) proposed a regression model approach to estimate the average causal effect of restricted mean survival times. Xie and Liu (2005) developed the adjusted Kaplan-Merier estimators of treatment-specific survival functions using inverse probability weighting. Stürmer et al. (2005) demonstrated that different strategies are comparable to adjust for confounding for estimating hazard ratios including inverse probability weighting, propensity score matching, and stratification for a given dataset. Zhang and Schaubel (2012) combined Chen and Tsiatis (2001)'s regression method and the inverse probability weighted Nelson-Aalen estimator resulting in a doubly robust estimator of the average causal effect of restricted mean survival times. Tchetgen and Robins (2012) also proposed inverse probability weighting and doubly robust methods for the causal hazard ratio.
The regression adjustment is susceptible to bias due to model misspecification; while inverse probability weighting methods are unstable when the propensity score is close to zero or one. Matching methods have been demonstrated to overcome these drawbacks and present many practical advantages Stuart;2010). Therefore, matching is widely used in observational survival studies 2007, 2008. However, there is no gold standard approach to analyze survival data after propensity score matching with replacement, and variance estimation is open to debate. Gayat et al. (2012) used the Cox proportional hazards model and the robust variance estimator to infer the causal hazard ratio accounting for the matched pairs. Therneau and Grambsch (2013) suggested fitting the Cox proportional hazards model stratifying by the matched pairs. However, Austin (2013) pointed out that this strategy may give biased results for the marginal hazard ratio due to non-collapsibility. On the other hand, ignoring the matched nature is inappropriate 2014). More importantly, in practice, the propensity scores are usually estimated through a parametric model in the first stage. Existing methods rarely incorporate the uncertainty prior to or in the matching stage and restrict inference to the matched sample.
In this article, we consider propensity score matching for causal inference in the context of survival analysis with two treatments. We focus on estimating the population average causal hazard ratio, which is arguably the parameter of most interest given the widespread use of the Cox model in practice. We note that paired matching may result in discarding data when the numbers of treated and control individuals are not equal. Following Abadie and Imbens (2006), we use propensity score matching as a tool for imputing unobserved potential outcomes, i.e. counter-factual outcomes, for all individuals. Following most empirical literature, we use a logistic regression model to estimate the propensity score function. Because the propensity score function is estimated prior to matching, it is necessary to account for the uncertainty due to parameter estimation. Due to the non-smooth nature of matching, our derivation is based on the technique developed by Andreou and Werker (2012), which offers a general approach for deriving the limiting distribution of statistics that involve estimated nuisance parameters. This technique has been successfully used in Abadie and Imbens (2016) for the matching estimators of the average causal effects in a continuous outcome context. We extend their results to the matching estimator in the survival outcome context. This extension is not trivial, because the survival outcome is often right censored. We utilize the martingale theory (Hall and Heyde;1980) of the counting process to establish the asymptotic distribution of the propensity score matching estimator of the causal hazard ratio.
Lack of smoothness also makes the conventional replication methods invalid for variance estimation for the propensity score matching estimator. The bootstrap method (Efron;1979) is a commonly used resampling approach that permits estimation of the sampling variability of estimated parameters. There is some controversy surrounding variance estimation when estimating treatment effects in propensity score matched samples. Abadie and Imbens (2008) demonstrated the failure of the bootstrap variance estimator for two-treatment matching with replacement and a fixed number of matches. The main reason is that the bootstrap sample does not preserve the distribution for the number of times that each individual is used as a match. Yet, Austin and Small (2014) demonstrated by simulation that the bootstrap has a good coverage property based on propensity score matching without replacement. Adusumilli (2017) recently proposed a bootstrap inference for the propensity score matching estimator for the continuous outcome, which improves the variance estimator based on the asymptotic distribution described by Abadie and Imbens (2016). Accordingly, we propose a new replication-based inference for the propensity score matching estimator in the survival outcome context. Simulation demonstrates that the replication-based inference outperforms Wald-type inference based on the asymptotic theory in finite samples.

Potential outcomes
We use the potential outcomes framework 1974). We assume that the individuals are random samples from a larger population. Let X ∈ R d be a vector of baseline covariates. Let W denote the treatment, which takes value in {0, 1}. Let T (ω) denote the potential time to a clinical event had the individual received treatment ω. Let C (ω) be the time to censoring had the individual received treatment ω. We assume non-informative censoring under which C (ω) ⊥ ⊥ T (ω) , where ⊥ ⊥ means "independent of". The full set of variables is F = (X, W, T (0) , T (1) , C (0) , C (1) ).

Marginal structural Cox proportional hazards model
For comparing two treatments, we consider the ideal case where we can observe F for each individual. We define the causal parameter as the hazard rate of failing at time t for a population of patients had they received treatment ω. We define U (ω) = min(T (ω) , C (ω) ) and ∆ (ω) = I(T (ω) ≤ C (ω) ) to be the time to a clinical event or censoring and the clinical event indicator, respectively, for an individual received treatment ω. We also define the counting process N (ω) (t) = I(U (ω) ≤ t, ∆ (ω) = 1) and the at-risk process Y (ω) (t) = I(U (ω) ≥ t).
Following Tchetgen and Robins (2012), we assume a causal Cox proportional hazards model.
Definition 1 (Causal Cox Proportional Hazards Model). The structural model for comparing treatments is where λ 0 (t) is the population hazard rate if all individuals had treatment ω = 0.
This is a causal model, because it compares the outcome under different treatments in the same group of individuals, i.e., the population of all individuals. The parameter β describes the relative hazard of having a clinical event if all individuals received treatment ω = 1 compared to if all individuals received treatment ω = 0.

Observed data and the no unmeasured confounding assumption
The observed data are summarized as where X i is a vector of baseline covariates, W i is the actual treatment, T i is the time to a clinical event, C i is the observed censoring time, U i is the time to a clinical event or censoring, ∆ i is the clinical event indicator, We now define the observed data counting process as N i (t) = I(U i ≤ t, ∆ i = 1) and the at-risk process as Y i (t) = I(U i ≥ t). Let e(X i ) = P (W i = 1|X i ) be the propensity score for individual i to receive treatment ω = 1.
We make the consistency assumption that links the observed data processes with the potential outcome processes.
In order to use the distribution of the observed data to estimate the parameters in Model (1), we require the assumptions of unconfoundedness (also known as no unmeasured confounding or treatment ignorability) and positivity (Robins;. Assumption 3 (Positivity). There exist positive constants c andc such that with probability 1, 0 < c < e(X i ) <c < 1.

Matching
If we were to observe all potential outcome processes, we can fit a Cox proportional hazards model to the data {N : ω = 0, 1; i = 1, . . . , n} to obtain an estimator for the causal log hazard ratio β 0 . However, the fundamental problem in the potential outcome framework is that we can observe for a particular individual i, {N (t)} with either ω = 0 or ω = 1 but never both. From this point of view, causal inference is inherently a missing data problem.
We use matching as a tool to impute the missing potential outcomes. We consider matching with replacement with the number of matches fixed at one, ignoring ties. To be precise, for individual i, the potential outcome process under W i is the observed outcome process (U i , ∆ i ) ; the potential outcome process under 1 − W i is not observed but can be imputed by the observed outcome process of the nearest individual with 1 − W i . Table 1: Imputed data structure: k i is the number of times individual i is used as a match.
We first consider the case when the propensity score is known and matching is done based on the true propensity score. Define the propensity score matching function m(ω, p) as the index of the individual with treatment level ω that is closest to p in terms of the propensity score values m(ω, p) = arg min Without loss of generality, we use the Euclidean distance to determine neighbors; the discussion applies to other distances 2006). Define the imputed potential outcome as or equivalently the imputed potential outcome process as for ω = 0, 1. Table 1 presents the full imputed dataset. Let k i (W i ) be the number of times that individual i is used as a match for the opposite treatment group. Then, the full imputed dataset can be equivalently represented by a weighted observed dataset, where each individual i is assigned with a weight 1 + k i (W i ), 1 is the weight for individual i, and k i (W i ) is the weight individual i can represent for individuals in the opposite treatment group.

Estimating equations
Define Λ 0 (t) = t 0 λ 0 (v)dv as the cumulative hazard function at time t. In what follows, we derive our estimators for β 0 and Λ 0 (t), t ≥ 0.
Based on the imputed dataset, the estimating functions for β 0 and Λ 0 (t), t ≥ 0 are where (3). We can also write (4) and (5) as Setting (6) equal to zero, we can obtain the estimator for dΛ 0 (t) for fixed β as Substituting (8) into (7), we obtain the estimating equation to solve for β. Let then β is the solution to the estimating equation which is the partial score equation for a proportional hazards model with W i as the covariate and weights 1 + k i (W i ) for each individual i. Therefore, β can be calculated by standard software.
Theorem 1. Under the regularity conditions presented in the supplementary material, with the known propensity score, and In observational data, the propensity score is often unknown and therefore has to be estimated from the observed data. Following most of the empirical literature, we assume that the propensity score follows a generalized linear model, e(X T i θ). We can then obtain a consistent estimateθ. The matching procedure can be carried out with the estimated propensity score e(X T iθ ). We now denote k i (ω) to be kθ ,i (ω) to reflect its dependence on the estimateθ.
To summarize, the algorithm for developing a matching estimator for β 0 when the propensity score is unknown is the following.
Step 1. Using the data {(X i , W i ) : i = 1, . . . , n}, fit a propensity score model, e.g. a logistic regression model e(X T i θ) and obtain an estimateθ.
Step 2. Based on the estimated propensity scores, create an imputed dataset as in Table 1, where the missing potential outcomes are imputed using the nearest neighbor imputation.
Step 3. Create weights for each individual i as 1 + kθ ,i (W i ), where kθ ,i (W i ) be the number of times that individual i is used as a match for the opposite treatment group based on the estimated propensity score estimates.
Step 4. Obtain estimator for β 0 by solving (10) using the standard software; e.g., the function "coxph" in R with the weighting argument.
We study the asymptotic properties of the estimator in the case when the matching variable is the estimated propensity score. The technique we will use is based on Andreou and Werker (2012), who proposed a unified approach to study the limiting distribution of statistics with estimated nuisance parameters. This technique has been applied by Abadie and Imbens (2016) to study the asymptotic distribution of the matching estimator based on estimated propensity score in the context of causal inference.
Theorem 2. Under the regularity conditions presented in the supplementary material, with the estimated propensity score, Theorem 2 shows that matching based on the estimated propensity score improves the efficiency of the matching estimator compared to matching based on the true propensity score if it is known. This phenomenon is in line with that in the setting with a continuous come .
Theorem 2 provides variance estimation based on approximation of the asymptotic variance formula. See the supplementary material. In the next section, we propose re-sampling variance estimation that has a better finite-sample performance. Abadie and Imbens (2008) demonstrated that the nonparametric bootstrap is invalid for the matching estimator based on matching with replacement and with a fixed number of matches. Otsu and Rai (2017) suggested resampling the matching estimator directly based on its martingale residual terms, which only works for matching based on the full vector of covariates. In order to reflect the uncertainty in the estimation of the propensity score, Adusumilli (2017) proposed a hybrid bootstrap procedure by resampling treatment variables and resampling the martingale residuals under both treatment conditions. This necessitates imputation of the unobserved martingale residuals under the opposite treatment. We define the martingale residuals as

Resampling based inference
where µ H (ω, p) is obtained by a non-parametric regression estimation of H i on e(X T i θ) among individuals with W i = ω. We propose a double-resampling procedure as follows.
Step 0. For individual i, find the secondary nearest neighbor matching function m(ω, X i ); if W i = ω, m(ω, X i ) = i, otherwise, m(ω, X i ) returns the index of the nearest neighbor in the opposite treatment group, matching based on the full X, rather than on the propensity score. The pair of potential residuals for individual i are Moreover, we need to estimate the matching function we use the following imputation strategy: create quintiles based on e(X T i θ), identify the quintile individual i falls, randomly sample one, say l, from that quintile with W j = ω, and let k i (ω) = k θ,l .
Step 1. Generate the bootstrap treatment, W * i , i = 1, . . . , n, from Step 2. Based on (W * i , X i ) n i=1 , re-fit the propensity score model and obtain the replicate θ * .
Step 3. Generate a sequence of independent and identically distributed random variables (u * 1 , . . . , u * n ) from u * i ∼ N (0, 1) or a two-point discrete distribution Step 4. Define the new bootstrap residuals as whose expectation over the probability distribution implied by Step 1, conditional on the original data, is where k θ,i (ω) is imputed at Step 0. Re-center the bootstrap residuals and compute the Step 5. Repeat Steps 1-4 for a large number B times and denote the bth replicate of S n (β 0 ) as S * (b) We construct the (1 − α)% confidence interval for β 0 as ∂S n ( β)/∂β −1 S * n(α/2) , S * n(1−α/2) . Remark. Here we provide some discussion to Step 0. Note that the secondary matching procedure matches on the full set of covariates, rather than on the propensity scores. Doing so preserves the conditional covariance in (16) between X and the error terms r 2i , given the propensity scores. This term reflects the improvement when using the estimated propensity score. See Adusumilli (2017). We impute k i (ω) by drawing a value from the empirical distribution of k i (W i ) in propensity score strata for W i = ω to re-create the distribution of k i (W i ).

Simulation
We conduct simulation studies to assess the finite sample performance of the proposed estimators relative to existing estimators: (i) the unadjusted estimatorβ nai , which is obtained by fitting the Cox proportional hazards model with the treatment status as the only covariate and therefore is used to demonstrate the degree of confounding bias; (ii) the inverse probability weighting estimatorβ ipw , where a weighted Cox proportional hazards model was fitted with the treatment status as the only covariate, weighted by the inverse of the probability of the actual treatment; (iii) the regression estimatorβ reg , where the Cox proportional hazards model was fitted with baseline covariates and treatment status as the covariates; (iv) the propensity score matching estimatorβ psm,0 based on the true propensity score; (v) the proposed propensity score matching estimatorβ psm based on estimated propensity score.
These comparative estimators have been widely used in practice, which motivates us to compare these estimators in our simulation study. For variance estimation ofβ nai ,β ipw , andβ reg , we consider the robust output from the standard software, denoted byV nai ,V ipw , andV reg . For inference based onβ psm , we construct Wald confidence interval based on the empirical variance estimator and two bootstrap percentile confidence intervals with the conventional bootstrap and the proposed double-resampling procedure.
We consider a sample size of n = 1000. Let X i = (X 1i , X 2i ) be a vector of baseline covariates, where we generate X ki from exponential distribution with mean 1/λ k with λ k = 1 for k = 1 and 2. We generate the treatment indicator W i ∼Ber(p i ), where p i = {1 + exp(θ 0 + θ T 1 X i )} −1 . We generate two potential outcomes T (0) i and T (1) i according to Algorithm S1 that are compatible with Model (1). The observed time to event is We also generate the time to censoring C i from a uniform distribution, so that 20% to 30% of the event time is censored.
We consider simulation scenarios formed by crossing the following two factors: 1) degree of confounding ranging from weak, medium and strong by considering θ = (θ 0 , θ 2) model specification of the propensity score model: Scenario (i) using a correct specification of the propensity score model; i.e., a logistic regression model with the linear predictor θ 0 + θ T 1 X i , and Scenario (ii) using a misspecification of the propensity score model; i.e., a logistic regression model with the linear predictor θ 0 + θ T 1 X 1/2 i . Table 2 shows the results when true β 0 = 0. The results for β 0 = −0.5 and β 0 = 0.5 are presented in the supplementary material. The unadjusted estimatorβ nai has severe bias and thus barely captures β 0 , which becomes worse for stronger separation of propensity scores between treatment groups. The inverse probability weighting estimatorβ ipw greatly reduces the bias; however, its variability is underestimated when two treatment groups are strongly separated in propensity scores and the propensity score model is correctly specified, thus the robust confidence interval from the standard software has some undercoverage issue in this case. The bias of the regression estimatorβ reg is always much larger than the propensity score based estimators. For all investigated simulation scenarios,β psm possesses the smallest bias and best coverage rate of 95% confidence intervals using the proposed asymptotic variance estimation and double-resampling variance estimation, among all estimators under comparison. On the other hand, using the standard software's robust method always overestimates the variances ofβ psm,0 andβ psm . Thus, our newly proposed variance estimation approaches is apparently beneficial for making reliable inference. Simulations with a sample size of n = 200 are included in the supplementary material. For the small sample cases, it is noticed that the double-resampling method recovers the 95% confidence level better than the asymptotic method.

An application
Non-small cell lung cancer (NSCLC) is the most commonly diagnosed lung cancer; typically, around half the NSCLC patients receiving chemotherapy will receive additional treatment in the post progression setting, i.e. second-line treatment setting, where "carboplatin + paclitaxel" and "erlotinib" are two historically commonly used treatments; see Cui et al. (2018). In this section, we use the IMS Health Oncology electronic medical record (EMR)  Note: "Var" is the variance of point estimates of β 0 across simulated datasets; "VE" is the average variance estimation for the point estimators over simulations, thus VE minus Var reflects the bias in estimated variance; "CR" is the empirical coverage rate of 95% confidence intervals. Four types of variance estimates forβ psm were compared: "software", output from the standard software; "asymp", the proposed asymptotic variance estimation; "naiveboot", the naive nonparametric bootstrap; "double-rsp", the proposed double-resampling method. data to conduct a comparative effectiveness analysis of the two treatments in the secondline setting with the proposed causal inference of hazard ratio based on propensity score matching estimator and other existing estimators mentioned in the previous section.
The EMR data is deidentified observational patient-level clinical data with demographic and baseline clinical characteristics collected from medium and large community-based oncology practices across 50 states of the USA. The dataset used contains retrospective cohort of 10,634 eligible patients at least 18 years old who received at least two lines of therapy, from 1 January 2007 to 31 December 2014; see Cui et al. (2018) for details. The time from the start date of second-line therapy to the death date was defined as overall survival. Patients alive at the end of the time period were censored at the end date of the dataset, and patients with the time between last visit and the end of dataset shorter than a sufficient time were censored at the date of last visit, where sufficient time is defined as twice the average time between visits in 3 month before the last visit. Missing data were classified into its own category for each categorical variable; see de Rooij (2018). Among the eligible patients in the dataset, 1241 patients were treated with carboplatin + paclitaxel, and 895 patients received single-agent erlotinib as second-line therapy.
The propensity scores were estimated using a logistic regression model with predictors: age at initiation of second-line therapy, gender, race, region, disease stage at initial diagnosis, Eastern Cooperative Oncology Group performance status score at initiation of second-line therapy, facility types of academic or community cancer center, year of index diagnosis, and days from index diagnosis to initiation of second-line therapy, see Cui et al. (2018) for details. To improve the overlapping in the covariate distribution for more credible matching, see Yang et al. (2016), we applied Crump et al. (2009)'s trimming method to drop patients with extremely high or low propensity scores, the remaining sample contains 1027 patients receiving carboplatin + paclitaxel and 872 patients receiving erlotinib. We proceed with the data of the remaining patients in the two cohorts for comparative effectiveness analysis on the population represented after trimming. Table 3 shows estimated log hazard ratioβ, point estimates and 95% confidence intervals for the hazard ratio of "carboplatin + paclitaxel" to "erlotinib" based on the unadjusted estimatorβ nai , the inverse probability weighting estimatorβ ipw , the regression estimatorβ reg with the robust variance estimation from the standard software for constructing Wald confidence intervals. For the proposed inference based on propensity score matching estimatorβ psm , we constructed both Wald confidence interval based on the empirical asymptotic variance estimator and bootstrap percentile confidence interval with the proposed double-resampling method. All adjusted methods gave higher hazard ratio estimates than the unadjusted naive method. The propensity score matching approach estimated that "carboplatin + paclitaxel" reduces the hazard by around 3.5% versus "erlotinib", with both the 95% confidence intervals based on the proposed asymptotic method and double-resampling method including 1 for hazard ratio, indicating insufficient evidence for a statistically significant difference at 0.05 level in the effectiveness comparison. Among the other adjusted approaches presented in Table 3, another 95% confidence interval covering 1 for hazard ratio is based on inverse probability weighting. The proposed asymptotic method and double-resampling method for propensity score matching estimator provided wider intervals than other adjusted approaches in Table 3, while the double-resampling method seemed to be more conservative than the asymptotic results in effectiveness comparison.

Discussion and future studies
We establish the statistical properties of the propensity score matching estimator of the marginal causal hazard ratio based on matching with replacement and a fixed number of matches. We also propose a double-resampling technique for variance estimation that takes into account the uncertainty due to propensity score estimation prior to matching. Several issues are worth discussing. First, although the literature has proposed flexible machine learning methods for estimating the propensity score function, our development relies on a parametric assumption, which enables to quantify the impact of propensity score estimation on the matching estimator. Second, there may not be sufficient overlap in terms of distribution of covariates between treatment groups. Crump et al. (2009) and Yang and Ding (2018) suggested trimming the propensity score beyond some thresholds. The criterion is based on minimizing the asymptotic variance of continuous outcome. It would be interesting to develop similar criteria to identify optimal overlap regions for the survival outcome. Third, Yang et al. (2016) has developed the generalized propensity score matching for multiple treatments for continuous outcomes. It is also of interest to extend our inference framework to more than two treatment comparison for survival outcomes. Fourth, in our future work, we will develop sensitivity analyses to no unmeasured confounding and methods to handle missingness and informative censoring.

Supplementary material
Supplementary material includes proof.

S1 Regularity conditions and Lemmas
In this section, we provide the regularity conditions and lemmas.
Assumption S1 (Regularity conditions). (i) X has a continuous distribution with a compact support and a continuous density function; The regularity conditions are standard in the literature and are often assumed for technical convenience.
Lemma S1. Suppose that (W 1 , X 1 ),...,(W n , X n ) are independent and identically distributed, where X is a scalar continuous variable with a bounded support. Suppose also that σ 2 H (ω, x) is uniformly bounded over the support for X. Then, for a given ω, in probability, as n → ∞.
S2 Proof of the asymptotic unbiasedness of n −1 S n (β 0 ) This section includes three parts that follow the similar logic of proof. The first and the second parts provide some results useful for later sections. The proof for the asymptotic unbiasedness of n −1 S n (β 0 ) is located in the third part.
For ω = 0, 1, define dM (ω) (t) = dN (ω) (t) − dΛ 0 (t) exp(β 0 ω)Y (ω) (t). From the standard theory for the counting process, dM (ω) (t) is a martingale process with respect to the population and its baseline hazard Λ 0 (t). Next we will prove that I( } is a martingale for the imputed pseudo-population which means that imputed pseudo-population has similar covariates distribution with the target population. First, we show that for ω = 0, 1, as n → ∞. We show (S1) for ω = 1. The proof for ω = 0 is similar and therefore omitted. We express (S1) for ω = 1 as where the second line follows by the consistent assumption, and Under Assumption S1, Abadie and Imbens (2006) showed that k i (1) δ is bounded almost surely for any δ > 0, and the discrepancy due to matching is |X m{1,e(X i )} −X i | = O p (n −1 ) for a scalar X. It follows that T 1 is consistent for zero. Moreover, assuming that E dM (1) (t) | p is Lipschiz continuous in p, T 2 is consistent for zero. Lastly, by the strong law of large numbers, T 3 is consistent for zero. Therefore, (S1) follows. Since

S3 Proof for Theorem 1
Taylor expansion of S n ( β) = 0 around β 0 leads to whereβ is on the line segment between β and β 0 . Assuming that n −1 ∂S n (β)/∂β T is invertible, we have Under suitable regularity conditions, n −1 ∂S n (β)/∂β T converges in probability to a positivedefinite matrix A(β 0 ). Then the reminder is to show the asymptotic distribution of n −1/2 S n (β 0 ).
Theorem S1. Suppose Assumptions 1-3 and Assumption S1 hold. Then, in distribution, as n → ∞, where Proof. We will show that n −1/2 S n (β 0 ) can be expressed as a sum of n independent and identically distributed random vectors plus a term that converges in probability to a zero vector. By some algebra, we obtain Therefore, continuing with (S3), we obtain where (S10) follows from (S8) . Under the standard regularity conditions in survival analysis, in probability, as n → ∞. Therefore, (S11) becomes where H i (ω) is defined in (11). Moreover, where the last line follows by the martingale property for the potential outcome process.

S5.1 An data-generating algorithm
We describe the algorithm for generating T (0) and T (1) that are congenial with Model (1) as follows.

S5.2 Illustration of propensity score distributions with weak, medium, and strong separation
This section demonstrates the mentioned simulation settings of weak, medium, strong separations in the propensity scores between the two treatment groups. Density curves for the true propensity scores of the two treatment groups are presented in dashed lines of W = 0 and in solid lines of W = 1 with weak, medium, strong separations respectively in Figure  S1.

S5.3 Additional simulation results
We provide additional simulation results. The simulation setting are the same as in §6 except that β 0 is set to be −0.5 and 0.5.  Figure S1: Plots of weak separation, medium separation, and strong separation in distributions of propensity score between treatment groups W = 0 and W = 1. Note: "Var" is the variance of point estimates of β 0 across simulated datasets; "VE" is the average variance estimation for the point estimators over simulations, thus VE minus Var reflects the bias in estimated variance; "CR" is the empirical coverage rate of 95% confidence intervals. Four types of variance estimates forβ psm were compared: "software", output from the standard software; "asymp", the proposed asymptotic variance estimation; "naiveboot", the naive nonparametric bootstrap; "double-rsp", the proposed double-resampling method.