“Lead Time impact on bullwhip effect in an Order-up-To-Policy single echelon supply chain with ARMA Demand ”

(1)

Master Thesis for Technology and Operations

Management Course, RUG University, Groningen

Student: Michele Pierangeli S3002470

Supervisor 1: Prof.Ir. G.J.C. Gaalman

Supervisor 2: Prof. Ir. R. Teunter

“Lead Time impact on bullwhip effect in an

Order-up-To-Policy single echelon supply

(2)

2

Abstract

The Bullwhip effect has been investigated under a myriad of perspective and assumptions in order to understand a way to hamper its tremendous effects on revenues and logistic systems: in this

paper, a single-echelon supply chain with order-up-to policy inventory system, and ARMA(n,n) demand forecast is enquired to understand the impact of lead time on bullwhip behaviour. Analytical and simulation analysis are performed in order to test the relationship of ARMA(n,n) order, poles, zeroes and the increment of Bullwhip effect as a function of lead time: by employing

Tsypkin theorem and Z transform, and a new quantification of the Bullwhip effect developed by Gaalman & Disney (2012), demand and orders variance of the system will be calculated through impulse response analysis. Finally, sufficient conditions for the increment of Bullwhip effect along

lead time are derived and proved. Managerial and practical implications are briefly discussed.

Introduction

Well begun, is half done.

Italian proverb

The Bullwhip effect phenomenon, or Forrester effect, is the amplification, propagation of variance upstream in the supply chain (Gaalman, 2006). More simply, the fluctuations of the quantities ordered from suppliers at each stage of the supply chain increase along the supply chain itself, from retailer to producers. This name was first utilized by P&G, and before reaching academic consensus was also known as whiplash or whipsaw effect (Lee et al. 1997). The phenomenon showed was first recognized and approached in the 50’s (see Forrester,1958, 1961, Magee & Boodman ,1967, etc.), but it is only later that simulation and games (Sterman 1989, Goodwin & Franklin, 1994) showed the bullwhip persistent arising, even in unstructured games. The

persistence of this effect can be dreadful: as noted by Metters (1997), the economic consequences of the bullwhip effect can reach the 30% of factory gate profits (Gaalman, 2005).

In this research, we will approach the Bullwhip effect utilizing an order-up-to policy to analyze its effects under a control engineering view: an impulse response analysis of the model will be performed utilizing the important results of Tsypkin Theorem which links the variance and the impulse response of a system (Tsypkin, 1958, 1964), and in this case will make calculations very simple. In fact, since the bullwhip effect compares the variances of the demand and the order quantities at a certain stage of (or along all) the supply chain, Tsypkin’s Theorem applies by expressing these variances through a linear expression of the system’s impulse response. In this research will exploit a different measure introduced by Gaalman and Disney (2012, working paper): instead of the classical ratio among order and demand variance, we will utilize their difference over the impulse response variance. This formulation, as shown later, will make it possible to exploit Tsypkin Theorem’s results, which allows us to express linearly the variance of the whole system.

(4)

4 position “-up-to” a pre-defined level (Gaalman, 2005 ). Anyway, due to a lead-time between placing an order and physically receiving the goods to stock, in case of uncertain demand, we have to make forecasts (Chen and Disney, 2003), in order to build a realistic model which can cope with variable demand.

To forecast the demand, Auto Regressive Moving Average models (generally known as ARMA models) will be utilized: ARMA models first definition has been developed by Box and Jenkins in 1970 (Box, Jenkins, 1970), and nowadays they are utilized for modelling in engineering and applied economics, due to their flexibility (Lucchetti, 2015). ARMA models describe a stochastic processes which is made of an autoregressive component (AR(p)) consisting of p parameters and a moving average part (MA(q)) of q parameters. Being more specific, this research is based on the work of Gaalman & Disney (2012, working paper), who introduced in the Amsterdam POM World

conference their preliminary analysis on the bullwhip effect with ARMA(2,2) modelling demand in an order-up-to policy. Thus, the main proposition we are going to enquire is the following,

proposed by Gaalman and Disney (2012, working paper).

Proposition 1: “If and only if the impulse response of CB(k) > 0, then CB(k) > 0 for every k>0 and is an increasing function of lead time k”

This proposition links the positivity of the impulse response of the system over (lead) time to the ratio of the order and demand quantities: by increasing lead time, a positive impulse response should incrementally (since the system is assumed stable, and will eventually head to a finite value) increase the bullwhip effect. In other words, the sum of positive values as indicated by Tsypkin’s theorem, will grow with the number of terms (or periods) in exam (lead time), if and only if all the values of the impulse response will be positive.

To ascertain this matter, we introduce proposition 2 and 3. The original expression of proposition 2 and 3 has slightly changed during the course of the research: they have been modified following what empirical results suggest. It is common knowledge in control engineering that a discrete system with poles (the solution of the polynomial denominator of the ARMA(n,n) polynomial expression) and zeroes (the solution of the polynomial numerator of the ARMA(n,n) polynomial expression) which are bigger than 1 (>1) are not stable. With negative poles and zeroes, the results are oscillatory. So we collocate our analysis region among 0 and 1: thanks to this, some orderings of the poles and zeroes of the ARMA(n,n) models in exam do always follow proposition 1 and some do never. Others have ambivalent behavior. A proper explanation will be shown in the following sections.

The condition is the following: given an ARMA(n,n) model (with poles and zeroes located on the real axe between 0 and 1), with n poles (or eigenvalues, 𝜆/_. ) and n zeroes (𝜆0.), which are ordered from the largest (with index i=1) to the smallest (with index i=n), we are interested in the

difference between , 𝜆/_. and 𝜆.0, for every i from 1 to n. If every difference is >0, bullwhip effect increases over lead time; if every difference is <0, will never be possible , whatever the value of 𝜆_./ and 𝜆0. (provided that ∀𝑖, 𝑗 𝑓𝑟𝑜𝑚 1 𝑡𝑜 𝑛: 𝜆𝑖𝜙≠ 𝜆𝑗𝜃):

So:

𝑷𝒓𝒐𝒑𝒐𝒔𝒊𝒕𝒊𝒐𝒏 𝟐:𝑔𝑖𝑣𝑒𝑛 𝑎 𝐴𝑅𝑀𝐴(𝑛, 𝑛) 𝑚𝑜𝑑𝑒𝑙 𝑤𝑖𝑡ℎ 𝑛 𝑧𝑒𝑟𝑜𝑠 𝑎𝑛𝑑 𝑛 𝑝𝑜𝑙𝑒𝑠 𝑚𝑜𝑑𝑒𝑙𝑙𝑖𝑛𝑔 𝑑𝑒𝑚𝑎𝑛𝑑 𝑎𝑠 𝑠ℎ𝑜𝑤𝑛 𝑏𝑦 𝐺𝑎𝑎𝑙𝑚𝑎𝑛 𝑎𝑛𝑑 𝐷𝑖𝑠𝑛𝑒𝑦 (2012) , 𝑔𝑖𝑣𝑒𝑛 𝑛𝑜 𝑝𝑜𝑙𝑒

(5)

5 𝑡ℎ𝑒 𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 𝑡𝑜 𝑜𝑏𝑡𝑎𝑖𝑛 𝑖𝑛𝑐𝑟𝑒𝑎𝑠𝑖𝑛𝑔 𝑏𝑢𝑙𝑙𝑤ℎ𝑖𝑝 𝑎𝑠 𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛 𝑜𝑓 𝑡ℎ𝑒 𝑙𝑒𝑎𝑑 𝑡𝑖𝑚𝑒 𝑖𝑠 𝒕𝒉𝒂𝒕 ∀𝒊 𝒇𝒓𝒐𝒎 𝟏 𝒕𝒐 𝒏: 𝜆/_. − 𝜆0_. _{> 𝟎.} 𝑷𝒓𝒐𝒑𝒐𝒔𝒊𝒕𝒊𝒐𝒏 𝟑: 𝑔𝑖𝑣𝑒𝑛 𝑎 𝐴𝑅𝑀𝐴(𝑛, 𝑛) 𝑚𝑜𝑑𝑒𝑙 𝑤𝑖𝑡ℎ 𝑛 𝑧𝑒𝑟𝑜𝑠 𝑎𝑛𝑑 𝑛 𝑝𝑜𝑙𝑒𝑠 𝑚𝑜𝑑𝑒𝑙𝑙𝑖𝑛𝑔 𝑑𝑒𝑚𝑎𝑛𝑑 𝑓𝑜𝑟𝑒𝑐𝑎𝑠𝑡𝑖𝑛𝑔 𝑎𝑠 𝑠ℎ𝑜𝑤𝑛 𝑏𝑦 𝐺𝑎𝑎𝑙𝑚𝑎𝑛 (2006), 𝑡ℎ𝑒 𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 𝑡𝑜 𝑛𝑒𝑣𝑒𝑟 𝑜𝑏𝑡𝑎𝑖𝑛 𝑖𝑛𝑐𝑟𝑒𝑎𝑠𝑖𝑛𝑔 𝑏𝑢𝑙𝑙𝑤ℎ𝑖𝑝 𝑜𝑣𝑒𝑟 𝑙𝑒𝑎𝑑 𝑡𝑖𝑚𝑒 𝑖𝑠 𝒕𝒉𝒂𝒕 ∀𝒊 𝒇𝒓𝒐𝒎 𝟏 𝒕𝒐 𝒏: 𝜆_𝑖𝜙− 𝜆_𝑖𝜃< 𝟎.

Proposition 2 and 3 define limit conditions for having a sure increasing bullwhip over lead time, or be sure not to: the disposition of poles and zeroes over the 0-1 line, say, the dimension of the poles and zeroes, that make their position on the line, has a direct link with bullwhip increment over time.

In the next section, Theoretical Background, we will describe the related theory, emphasizing the explanation of the order-up-to level inventory policy, ARMA models, and control theory approach. Methodology will provide the instruments of the analysis, which will be explained in the dedicated section. Finally, limitations and conclusions will be drawn.

Theoretical Background

Order-up-to policy:

Honesty is the best policy

Alison Berry Wilkinson, Berry Wilkinson Law Group, Berkeley

This subsection will explain the inventory policy used in the research: an inventory policy defines a range of rules to organize an inventory, which is the stock of any item or resource used in an organization. An inventory system is the set of policies and controls that monitor levels of inventory and determine what levels should be maintained, when stock should be replenished, and how large orders should be (Sianesi, 2012). We will utilize the Order-Up-To-Level policy, often called as the base stock policy, a very appealing policy in case of very high demand, under both a practical and theoretical perspective: it is simple, close to optimal and reflects to a great extent real world practices (Larsen and Thorstenson, 2008). As stated by Disney and Towill (2003), ““ Baganha and Cohen (1998) postulate a puzzling idea: “inventory management policies can have a destabilising effect by increasing the volatility of demand as it passes up to the chain”, whereas “one of the principal reasons used to justify investment in inventories is its role as a buffer to absorb demand variability”. In other words, inventories should have a stabilizing effect on material flow patterns. So how is it that market variability is amplified rather than dampened?””. They end up identifying the reason lying beneath this stunning contradiction as poorly designed inventories. Gaalman (2005) explains Order-up-to policy (OUT from here after) very clearly: ”Periodically, we review our inventory position and place an “order” to bring the inventory position “-up-to” a defined level.” In every period it happens the following: inventory level is reviewed and ordering decision is made at the beginning of the period, then the customers order is received and demand is realized and fulfilled at the end of the period. However, it takes (at least) one period to receive the order placed. Unsatisfied demand of a period is fully backordered. Two costs are

(6)

6 and s cost parameters which are constant over time (Chen and Disney, 2003). In this paper, we assume only a linear ordering cost. The objective is to minimize the long-run average total cost per period.

This policy is often set by the company to coordinate orders for multiple items from the same supplier, where setup costs may be reasonably ignored (Chen, Disney 2003).

In formulas, the one below is the inventory balance equation for OUT policy, made by the three discrete quantities (Gaalman, 2005):

i

_t+1+k

=i

_t+k

+o

_t

−d

_t+1+k

where i(t+k) is the inventory at period t+k (k periods after the current moment) , o the order decision at current moment t that will arrive at the beginning of period [t+1+k]) and dt+1+k is the demand during period [t+1+k]. When k=0 the system shows no physical ordering delay (lead time), even if effect of the sequence of events (ordering at the beginning of a period, receiving during - at least – the following period) within the ordering process brings a delay of one period.

Specifically, lead time refers to “the time span between the moment an order is placed by the buyer and the moment when ordered goods are delivered. It consists of the order processing time, production time and transportation time and is sometimes as long as several months”. (Ju et al., 2013) The higher the order up-to level, the better the service, which, when inventory

obsolescence is not an issue, is generally quite high (the service level is defined as the probability that demand is fulfilled by the on-hand inventory, Duc et al. 2008a). The key factors that

determine the amount of inventory needed, are the length of the replenishment lead time, the desired service level (fill rate or in-stock probability) and the demand uncertainty (Sianesi, 2012).

Minner (2003) reviews a variety of different investigations on order-up-to policy, such as capacitated order up-to-policy, emergency order-up-to policy, order-up-to policy in single and multiple echelon supply chains, with deterministic or stochastic lead time, and so on: this

witnesses the variety of academic interest in this simple though powerful inventory policy. Duc et al. (2008a), Luong et Phien (2006) and J. Dejonckheere et al. (2003), all utilize the order-up-to policy to investigate the bullwhip effect, either singularly, either comparing it to other inventory policies. In the next section we will explain how ARMA models are part of our research, and how they link to the order-up-to policy, ending up with a short review of the current literature which binds order-up-to policy, ARMA models and the bullwhip effect.

ARMA Models

“If you do not think about the future, you cannot have one” John Galsworthy, Nobel prize for Literature 1932, Swan Song (1928) Pt. II, Ch. 6

(7)

7 MA terms, was appropriately specified (Makridakis and Hibbon, 1995). This means that any series xt can be modelled as a combination of past xt values and/or random errors at time t+i.

Moving Average could be defined as a weighted linear combination of white noise variables. White noise is the most simple stochastic process to be imagined: in fact, it is a process with (at least) first and second statistical moment; these moments are constant through time (so the white noise is stationary) but do not give any memory to the process (Cryer, Chan, 2008). Moving

average processes are represented through this equation: 𝑦_b = 𝜇 + f(𝜃_.𝜂_bh.)

i

.jk

where 𝜇 𝑖𝑠 𝑎 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 (𝑤ℎ𝑖𝑐ℎ 𝑐𝑜𝑢𝑙𝑑 𝑎𝑙𝑠𝑜 𝑏𝑒 0);

Autoregressive models are, as their name suggests, regression on themselves (Cryer, Chan, 2008). Specifically, a p-th order autoregressive model AR(p) satisfies the equation:

𝑦_b = f(𝜙_.𝑦_bh.) q .jr + 𝜂_b . Together they give ∶ 𝑦_b = (𝜇) + f(𝜙_.𝑦_bh.) q .jr + f(𝜃_.𝜂_bh.) i .jk 𝐴𝑅𝑀𝐴(3,2) => 𝑦_b= 𝜙_r𝑦_bhr+ 𝜙_|𝑦_bh|+ 𝜙_}𝑦_bh}+ 𝜃_r𝜂_bhr + 𝜃_|𝜂_bh| 𝑤ℎ𝑒𝑟𝑒 𝜂bh. 𝑖𝑠 𝑡ℎ𝑒 𝑤ℎ𝑖𝑡𝑒 𝑛𝑜𝑖𝑠𝑒 𝑎𝑡 𝑡𝑖𝑚𝑒 𝑡 − 𝑖

In the methodology section, we will show how to link ARMA models, Z transform and the order-up-to policy. The next section briefly reviews what Z transform consists of. The final section of Theoretical Background will link all the argument here presented, through highlighting the very purpose of this research: analyze the bullwhip effect.

Z transform

Pantha Rei (everything changes)

Eraclito The the Z-transform was already known to Laplace, and basically consisted in the discrete

formulation of the Laplace Transformation (Bolzern, 2008). It was re-introduced in 1947 by W. Hurewicz et al., and gives a tractable way to solve linear, constant-coefficient difference

equations. It can be defined as two sided (ranging from -¥ to +¥) or one sided (ranging from 0 to ¥), as many integral formulation (Bolzern, 2008). The one sided expression of the z transform is the following: ∑€ (𝑓_•(𝑧h•

•jk )). It generally allows researchers to avoid the solution of difficult integro-differential equations, by dividing the time-shifted variable by the z operator as many times as much are the difference in time periods among the earliest and the latest variables (eg.: x(t)=3x(t+3) -> x(t)=3*x(t)*z^(-3)).

The Bullwhip effect

“The question that naturally arises is, if variable production is indeed costly, why does a homo economicus choose to bullwhip?”

(8)

8 In this section the Bullwhip effect will be explained. A good starting point to orientate among the huge literature around Bullwhip effect is to start from a good review: the latest up-to-date “The bullwhip effect: Progress, trends and directions” by Whang and Disney (2016), who explain how Bullwhip research is divided into the operational and the behavioral stream, into empirical or analytical research. Our focus, is evidently on the operational, analytical stream. A myriad of paper have been written in order to understand what actually causes Bullwhip effect: however, in their seminal paper of 1997, Lee, Padmanabhan and Whang brought the bullwhip effect to academic attention and proposed 4 operational causes of the bullwhip effect: 1) Demand nature and forecasting methods 2) Price, 3) (rational) shortage gaming and 4) order (& batching) policies. On the other hand, it is often said that problems in supply chains are “80% people centered and 20% technology centered (Andrasky, 1994) and that’s why another discerning claim should be made: by taking only a reductionist (operational) or a psychological (behavioral) perspective, we will never be able to completely explain the supply chain variance.

Regarding the harmful situations the bullwhip may bring the firm to, Carlsson and Fullér (2000) have further summarised the negative impacts as follows:

• Excessive inventory investments throughout the supply chain to cope with the increased demand variability

• Reduced customer service due to the inertia of the production/distribution system • Lost revenues due to shortages

• Reduced productivity of capital investment • Increased investment in capacity

• Inefficient use of transport capacity

• Production schedules missed more frequently (Chen, Disney 2003).

Therefore, how can we understand the bullwhip effect? For what concerns bullwhip modeling, Whang & Disney (2016) advise to consider 6 elements: demand nature, time delay, forecasting policy, ordering policy and information sharing mechanisms.The related literature generally assumed single (Gaalman, 2005 and Gaalman and Disney, 2006) or multiple echelon supply chain (Wang and Jia 2013 or Costantino et al. 2013), arbitrary (Gaalman and Disney, 2006) or sthocastic lead time (Michna et al. 2013, or Duc et al. 2008), stationary (Duc et al. 2007, Luong and Phien, 2006) or non stationary demand (Graves, 1999) several forecasting methods (Zhang, 2004 for instance studied the impact of diverse forecasting techniques on the bullwhip effect) and different inventory management policies (Mòran and Barras, 2006), generally tailored on the typology, fill rate, number and variety of stocked items.

(9)

9 effect.

Duc et al. (2008a) who investigates ARMA(1,1) model and the impact of lead time on the Bullwhip effect, utilized a minimum square error estimation technique for demand forecasting, and a base stock policy, to get the analytical expression of the order and demand variances and found out an irregular behavior of the bullwhip effect respect to lead time, when varying phi and theta

coefficients. The base stock policy is a special case of the order-up-to inventory policy: an order is placed at the beginning of each period in order to increase the inventory level up to a

predetermined level. Order-up-to level (which they call St, or also service level) is optimal in terms of total cost for inventory systems where there is no fixed ordering cost, while holding and

shortage costs are proportional to the volume of on-hand inventory and shortage, respectively (Nahmias, 1997; Zipkin, 2000). They find out that the bullwhip effect occurs only when the

autoregressive coefficient of the demand process is larger than the moving average parameter. To be more precise, the bullwhip effect does not always increase as lead time increases: in fact, it behaves differently depending on the lead time’s number of periods being odd or even.

Luong et Phien (2006) investigated a simple supply chain with one retailer and one supplier in which the retailer employs base stock inventory policy (which we explained above). They further develop bullwhip measure for the second-order autoregressive demand process and then extend the findings to the general autoregressive model. AR(2) process shows how bullwhip effect is an increasing function for increasing lead time for positive 𝜙r 𝑎𝑛𝑑 𝜙|, while if the second

autoregressive coefficient is negative, the bullwhip measure is oscillatory. However, they highlight how, for non-positive autoregressive coefficients, an increasing lead time is not always causing increasing (or the arise of) bullwhip effect.

Methodology

The best and safest method of philosophizing seems to be first to inquire diligently into the properties of things, and establishing those properties by experiments, and then to proceed more slowly to hypotheses for the explanation of them.

Isaac Newton After having modelled demand through ARMA models, an Impulse response analysis is required in order to apply Tsypkin theorem: to do this, it is convenient to express a time dependent

relationship through Z transform.

From ARMA expression to Z transformation.

Hereby, the decisive step of transforming an inventory model into a control engineering investigable system is performed.

(10)

10 •1 − 𝜙_r𝐵 … − 𝜙_q𝐵q_„𝑦

b…r = (1 − 𝜃r𝐵−. . . −𝜃i𝐵i)𝜂b…r

The z-transformation will allow us to express the following through a transfer function which will be deployed into a pole-zero representation: we will show how to express the impulse response as a function of pole and zeroes of the transfer function.

We assume that the demand can be described by the sum of a constant mean term and a normally distributed random variable it is not only our assumption but frequently made in literature (e.g.: Luong, Phien (2006)).

Therefore the demand d (being stationary) is modelled through a constant term and a stochastic process, zt+1:

Note a fundamental property is: .

We are interested in the demand variance, or in other words the fluctuation around a constant mean (since we consider demand as a stationary stochastic process) so we may ignore the constant term of the demand expression, and will model only the stochastic process. The so obtained expression will be transformed through Z transform and an impulse response analysis will be performed.

Box and Jenkins (1994) considered the operator z as the inverse of the shift operator B. Thus, starting from the previous (re-arranged):

𝑧_b…r= (1 − 𝜃r𝐵−. . . −𝜃i𝐵i) •1 − 𝜙_r𝐵 … − 𝜙_q𝐵q_„𝜂b…r So, applying the operator z:

𝑧_b…r =(1 − 𝜃r𝑧hr−. . . −𝜃i𝑧hi) •1 − 𝜙_r𝑧 … − 𝜙_q𝐵𝑧hq_„ 𝜂b…r

we obtain a Z-transformed version of our ARMA demand. Without loss of generality, we indicate as m the maximum among p and q (m=max(p,q)) and the subsequent step is to multiply both numerator and denominator expressed in a compact form by zm_{in order to get a polynomial} expression and finally a pole/zero representation.

The impulse response of can be written in (unilateral, single sided ) Z-transform (

(11)

11 where the eigenvalues corresponds with the poles in the Z-plane and the eigenvalues with the zeroes.

A usual simplification is using the Partial Fraction Extension Property, which holds for nominators with a degree less than m. Because of this we rewrite the transfer function to

and by this the nominator has the degree . Note .

And thus finally, by decomposing the two polynomial system by their roots, to: H(z) = ∏‹‰Œ•(‡hˆ‰Š) ∏‹_‰Œ•(‡hˆ_‰Ž)

, where 𝜆0_._{are called the zeroes of the system, and 𝜆} .

/_{are called the poles of the system. When}

H(z) is a rational function, i.e., a ration of polynomials in z, we call: 1. The roots of the polynomial numerator as the zeroes of H(z), and 2. The roots of the polynomial denominator as the poles of H(z).

Anyway, no poles of H(z) can occur within the region of convergence since the z-transform does

not converge at a pole: instead, the region of convergence is delimited by poles. (Bolzern, 2008).

𝐻(𝑧) − 1 = f • 𝑐. (𝑧 − 𝜆_./‘ ’ .jr with So we have:

From well-known table we have

which holds for

(12)

12

So and

for

Now

Or abbreviated expressed as .

Note the ( ) impulse response can be written as

or

In case of repeated (common) eigenvalues, the calculation of coefficients becomes:

𝐻(𝑧) − 1 = ∏ •𝜆/− 𝜆.0„ − ∏ •𝜆/− 𝜆. /_„ “ rjr “ .jr ∏ •𝜆/_{− 𝜆} . /_„ “ rjr = f 𝑟.(𝑛) •𝑧 − 𝜆_./„. “ .jr with 𝑟.(𝑛), the ith common pole, which is calculated by:

𝑟_.(𝑛) = 1 𝑛 − 𝑖!• 𝑑 𝑑𝑧•𝑧 − 𝜆/„ “ ∗ (𝐻(𝑧) − 1)— |_‡jˆ_Ž And finally, retransforming in time domain:

𝑝(𝑘 + 1) = f[𝑟_.(𝑛) ∗

“

.jr

𝜆/›h.…r_. ∗ 𝑘!

(𝑘 − 𝑖 + 1)!]

So it becomes clear that positive impulse responses depend on the values of the poles and zeroes of ARMA process.

The analysis will prosecute with the re-transformation to the time domain Z-1_{(H(z)) and the link to} Tsypkin theorem and the variances of the system, in order to investigate the Bullwhip behavior.

Impulse response analysis and Tsypkin Theorem

“Our impulses are too strong for our judgement sometimes”

Thomas Hardy, Tess of the D’Urbervilles

(13)

13 In this section we will deal with mathematical implant of this research. Gaalman and Disney (2012, working paper) express in formulas the following alternative Bullwhip measure:

where ∑oo is the order variance, ∑dd is the demand variance and ∑ is the impulse frequence response variance (of the white noise).

Conceptually, after modeling the demand and implementing the OUT (say, Order-up-to) policy, we will test our model in order to understand if lead time without expicitely calculating the variances

Σ

_dd

, Σ

_oo

, Σ

𝜂𝜂..

In fact, Tyspkin showed ( in the 1940’s – 1950’s) that the sum of squared impulse of the system response is equal to the system variance divided by the input variance (white noise with squared impulse response =1 at time t=0 ). Towill and Disney (2003), show how to apply this theorem to supply chain modelling by adapting Tsypkin results to supply chain modeling: in this way, they claim that we can express Bullwhip measure through a function of the impulse response of the order variance over the impulse response of the demand variance.

For the bullwhip effect criterion we need to know and . Respectively the sum of squared

impulse of the demand as the system output and the sum of squared

impulse of the orders as system output , both with the white noise as input.

The demand impulse response can be derived directly from ARMA process

Note a fundamental property is:

The impulse response of ARMA(n,n) demand can be then calculated, but becomes unconfortable and resource-expensive to enquire: therefore, an easier approach is needed. That’s where Z transform comes at hand. In fact an expression like the one below, becomes extremely complicated to express in terms of poles and zeroes.

Starting from p0 as p(0) and so on:

yt = - ∑A(k)y(t-k) + ∑B(k)e(t-k) ,

which gives a transfer function of:

H(z) = ∑B(k)z

-k

/(1+∑A(k)z

-k

)

This representation takes advantage of the ARMA finite complexity (Wei, 1974): 1)! The unknown input et causes the random behavior;

2)! The averaging section (MA) makes a finite response filter (FIR) 3)! The recursive section (AR) makes an infinite response filter.

Geraard Gaalman, my supervisor, (together with Disney, 2012), utilized space state method to calculate impulse frequence response to verify the arising of bullwhip effect after modelling it through space state or impulse frequence response. In this thesis, his alternative measure (developed by Gaalman & Disney) of the Bullwhip analysis will be employed:

instead of the classic :

where ∑oo is the order variance and ∑dd is the demand variance.

Thanks to technique such as Jury’s inner methods, the space state representation will lead to the analysis of the eigenvalues of our model, which at the present state of knowledge do not fully tell if a lead time increasing will definitely give birth to the Bullwhip effect.

Regarding the model representation, Stpehen Disney,in his lecture slides of the course Supply Chain

Modelling @ Cardiff Business School (2014) proposes a representation through block diagrams of the order-up-to policy supply chain, utilizing z transform:

Analytical modelling steps: a proposed timeline.

In this research, mastering the analytical tool is the hardest part: stochastic modelling, ARMA models, Box-Jenkins methodology, Jury’s inners methods, state space analysis, impulse frequence response analysis and pole-zeros analysis, Nyquist methods, and so on.

(14)

-14

Moreover for we have

However for large m this impulse response expression becomes rather demanding to analyze. Gaalman & Disney (2012, working paper), derived the following expression, which combine Tsypkin theorem and to the linear measure of Bullwhip effect presented by Gaalman & Disney (2006) and show its equivalence to the classic measure:

,

.

Instead of the ratio we will consider as a bullwhip criterion the variance difference: .

If we have a bullwhip.

The can be written as function of :

. Therefore applying Tsypkin results:

or, since we

expressed lead time L as L=k+1: .

For the bullwhip effect criterion we need to know and , respectively the sum of squared

impulse of the demand as the system output and the sum of squared impulse of the orders as system output. Both with the white noise as input.

By proving that the impulse response of demand is always positive, say, by proving proposition 1,

for every k>0, it is straightforward that the previous expression is

positively growing as lead time (L) grows .

Minding what we just defined, for the Impulse Response method (Bolzern, 2008) we assign at the input function u(t) the values of a impulse function. The impulse function is defined as follows:

1 k m> + 1 1 2 2 1 1 k k k m k m m k m p =

f

p _- +

f

p _- + +

f

_- p _{- -} +

f

p -1( var ) oo 1 dd BI output input iance =æ_çS ö_÷>

S

è ø

ˆ, ˆ

2( var ) 1

o o o o BI input ouput forecast error iance

(15)

15 Imp(k) = •1 𝑝𝑒𝑟 𝑘 = 0_{0 𝑝𝑒𝑟 𝑘 ≠ 0 , k>0}

(16)

16

Analysis

In this section proposition 2 and 3 will be enquired. Prop 1 is proved for all cases (a necessary and sufficient condition to have for all lead time an increasing bullwhip effect). Prop 2 and 3 are sufficient conditions (thus not necessary). All cases are based on real poles and real zeroes between 0 and 1. This limitation will be discussed in the first subsection. In fact, a sufficient and

necessary formulation for Proposition 2 and 3 has not been found yet, neither by Gaalman and

Disney, neither by me, but empirical results actually confirm proposition 2 and 3.

The first subsection will make us familiar with the mathematical expression of the analysis

performed: we will show a change of approach and introduce “proposition 0”, a necessary starting point for our analysis. In the second subsection, ARMA(2,2) (only in case of common poles, since for distinct poles it has already been analyzed by Gaalman & Disney (2009)) and ARMA(3,3) in case of common poles will be fully analyzed. The third subsection will extend the results in case of ARMA(n,n), starting from ARMA(3,3) with distinct poles and showing the practical utility of the change of approach (presented in the first subsection) for high order ARMA models.

Limitations and Conclusion will terminate the paper.

Terminology and problem presentation

In the first subsection, we will introduce the terminology and describe how a change in approach towards the z-transformed ARMA(n,n) models will allow a series of statements (analytically or empirically supported) regarding the Bullwhip increment along lead time.

The innovative change of approach in my research lies in a shift of perspective regarding a

complex ARMA(n,n) model: stimulated by the propositions suggested in Gaalman & Disney (2012, working paper), I started analyzing the different pole-zeroes disposition along the axis 0-1. The choice of taking values among 0 and 1 has been made for stability reasons: first, because for unstable systems the variances do not exist (Crier, Chen, 2008). Thus, case of real eigenvalues

(poles and zeroes) the condition is simply . Regarding the second and the third

condition, let’s consider the root locus in exam: for negative or complex values of the poles, we may have alternating (seasonal and cycling) oscillating values of the demands , with which the OUT policy don’t cope well.

Gaalman & Disney (2012, working paper) original propositions dealt with how many poles being bigger than how many zeroes. I ended up conceptually breaking down the various dispositions in 0-X or X-0 blocks, with X symbolizing a pole, and 0 symbolizing (obviously) a zero. This simple change in perspective will do a long way, as shown later.

My reasoning started from Proposition 1: in order to have increasing Bullwhip effect, the impulse response CB(k) needs to be always positive. Therefore, by

re-transforming the ARMA(n,n) model into the time domain,

1 i , i 1 f q l l

- < <

(17)

17 we can check for its positivity at t=0 and at t=infinite. Having chosen to express the discrete

formulation of our impulse response utilizing k=t+1, we will consider the values of the impulse response function in k=0+1=1 and in k=1+infinite=infinite. Regarding k=infinite, being it a linear combination of exponential lower than one, the retransformed system will converge to 0+_{or 0}-_, being the dominant pole (say, the biggest pole) multiplied by a positive or negative coefficient. Regarding k=1 (or p(1)), it is straightforward to prove that the value of the function in k=1 (t=0) is the sum of zeroes subtracted to the sum of poles (proof below, in this page, definition 1.0). This condition suggested that (as stated in proposition 2) to obtain a positive impulse response, a necessary condition would be that each pole had to be bigger than each zero. Breaking down the pole-zero disposition brought to the proof of Proposition 0.

Before discussing these statements, let’s explain here the terminology adopted from now on. The Z transform

- The nth_{largest pole (which is the smallest) and n}th_{largest zero (which is the smallest) are} referred to as 𝜆“/𝑎𝑛𝑑 𝜆“0. Therefore the (first) biggest pole is indicated as 𝜆r/and the (first) biggest zero is indicated as 𝜆r0. In case of common poles we will use the notation # to graphicly indicate it, and

Figure 2 - 0-X blocks and pole-zeroes numeration Figure 3 - X-0 blocks and pole-zeroes numeration

- 𝑟. = 𝑡ℎ𝑒 𝑐𝑜𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡 𝑜𝑓 𝑖bž 𝑏𝑖𝑔𝑔𝑒𝑠𝑡 𝑝𝑜𝑙𝑒, 𝑖𝑛 𝑐𝑎𝑠𝑒 𝑜𝑓 𝑛𝑜 𝑐𝑜𝑚𝑚𝑜𝑛 𝑝𝑜𝑙𝑒𝑠 - 𝑟.,•(𝑛) = 𝑡ℎ𝑒 𝑐𝑜𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡 𝑜𝑓 𝑖bž_{𝑐𝑜𝑚𝑚𝑜𝑛 𝑝𝑜𝑙𝑒 𝑜𝑓 𝑡ℎ𝑒 𝑗}bž_{𝑏𝑖𝑔𝑔𝑒𝑠𝑡 𝑝𝑜𝑙𝑒 𝑖𝑛 𝑎𝑛 𝐴𝑅𝑀𝐴(𝑛, 𝑛) 𝑚𝑜𝑑𝑒𝑙} 𝑤𝑖𝑡ℎ 𝑠𝑜𝑚𝑒 𝑐𝑜𝑚𝑚𝑜𝑛 𝑝𝑜𝑙𝑒𝑠 - 𝑟.(𝑛) = 𝑡ℎ𝑒 𝑐𝑜𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡 𝑜𝑓 𝑖bž_{𝑐𝑜𝑚𝑚𝑜𝑛 𝑝𝑜𝑙𝑒 𝑖𝑛 𝑎𝑛 𝐴𝑅𝑀𝐴(𝑛, 𝑛) 𝑚𝑜𝑑𝑒𝑙 𝑤𝑖𝑡ℎ 𝑎𝑙𝑙 𝑛 𝑐𝑜𝑚𝑚𝑜𝑛 𝑝𝑜𝑙𝑒𝑠}

The first condition is to have a positive p(1): as proved by Gaalman & Disney (2012, working paper), the sums of the poles and zeroes in case of no common poles, or all common poles, is equal to the value of p(1). In formulas, ∑“ 𝜆/_.

.jr − ∑“.jr𝜆.0 > 0 means p(1)>0, since working out the characteristic polynomials shows that:

, and , .

From the impulse response of the demand process we have , knowing that

(18)

18

( we get: ,and .

Thus definition 1.0 : p(1) = , which also holds for common poles or

eventually also for common zeroes, and for some common poles. With this in mind, we proceed to Proposition 0.

Proposition 0

Proposition 0 can be divided into 2 symmetrical statements, Proposition 0.a for positive p(1) and

Proposition 0.b for negative p(1).

Proposition 0.a

Proposition 0.a helps us identifying on the spot if a pole-zeroes disposition can ever be positive, simply by inspection. Let’s start identifying the 0-X or X-0 blocks for every i from 1 to n

((𝜆_./, 𝜆_.0) ∀𝑖 = 1: 𝑛). For instance:

Figure 4 - Proposition 0.a

Proposition 0.a : “If the difference between the biggest pole (𝜆_r/ indicated as X on the graph) and the biggest zero (𝜆0• indicated as 0 on the graph), regardless of the other poles and zeroes disposition, is greater than “hr_“ , (𝜆_r/− 𝜆_r0 _>“hr

“ ), we may eventually have increasing bullwhip along lead time, and we can’t affirm there surely won’t be one, since the p(1) will start with a positive value”.

Proof 0.a:

Let’s take 𝜀, 𝜂 > 0 and as small as possible (lim(𝜀, 𝜂) → 0…_{), with 𝜀 > 𝜂 > 0.} By considering 𝜆_r/ − 𝜆_r0 ₌“hr

“ + 𝜀, we take:

𝜆_r/ = 1 − 𝜂 and 𝜆.0 = (1 − 𝜂) − §“hr_“ + 𝜀¨ = r_“− 𝜂 − 𝜀 . Therefore, we want to prove that:

(19)

19 𝑝(1) = f 𝜆_./ “ .jr − f 𝜆_.0 “ .jr = (𝜆_./− 𝜆0_._{) + f(𝜆} .©“ / “ .j| − 𝜆_.0_{) > 0} The limit case consists in taking ª∑ (𝜆“ /_.

.j| − 𝜆.0„| as big as possible, and of course, negative. If the biggest (negative) sum of the other poles and zeroes won’t make p(1) fall below zero, than all the other pole-zeroes disposition won’t either.

The most negative ∑ ( 𝜆“ /_.

.j| − 𝜆.0) can be achieved by placing the poles all in 𝜆./ = 0 + 𝜂 (smallest poles allowed, since if the poles are equal to 0 there’s no meaning in calculating the impulse response) and the zeroes as big as possible 𝜆0‰ =

r

“− 𝜂 − 𝜀 (which coincide with the zero 𝜆0•).

Figure 5 - Proposition 0.a

Therefore we have: ∑ ( 𝜆“ _./ .j| − 𝜆.0) = (𝑛 − 1) ∗ (𝜆/| − 𝜆0|) = (𝑛 − 1) ∗ (𝜂 − (_“r− 𝜂 − 𝜀) = (𝑛 − 1) ∗ [2𝜂 + 𝜀 −_“r]. Thus: 𝑝(1) = f 𝜆_./ “ .jr − f 𝜆_.0 “ .jr = (𝜆_r/− 𝜆_r0_{) + f(𝜆} . / “ .j| − 𝜆0_._{) =}𝑛 − 1 𝑛 + 𝜀 + (𝑛 − 1) ∗ [2𝜂 + 𝜀 − 1 𝑛] After calculations: 𝑝(1) = (𝑛 − 1) ∗ 2𝜂 + 𝑛 ∗ 𝜀 > 0, (𝑏𝑒𝑖𝑛𝑔 𝜂, 𝜀 > 0), and so Proposition 0.a is proved.

Proposition 0.b

Proposition 0.b : “if the difference between the smallest zero (𝜆0« indicated as 0 on the graph) and the smallest pole (𝜆_“/ indicated as X on the graph) is greater (in absolute value) than ¬“hr_“ ¬, and it’s negative, we can NEVER have a bullwhip effect, because the p(1) will start with a negative value”.

(Remark: when 𝜆/_. = 0 we never have increasing bullwhip effect as function of the lead time). Proof 0.b :

(20)

20 The placement of the other poles and zeroes, is therefore limited: in order to maximize p(1) in this situation, the other n-1 poles should be put at maximum in 𝜆/_“ = 1 − 𝜂 , with 𝜂 > 0, since in 1 they would cause instability. The other zeroes can’t be placed in a point smaller than 𝜆0« position, since ,by doing so, 𝜆0«won’t be anymore the smallest pole (becoming 𝜆“hr0 ), invalidating the initial proposition. Therefore, the other zeroes can be placed at most in the same position as 𝜆0“.

Figure 6 - Proposition 0.b

We can express p(1) as:

f 𝜆_./ “ .jr − f 𝜆_.0 “ .jr = (𝜆/_“ − 𝜆0_“_{) + f(𝜆} . / “hr .jr − 𝜆_.0₎ But ∑“hr(𝜆_./

.jr − 𝜆.0) is just equal to (𝑛 − 1) ∗ (𝜆.©“/ − 𝜆0.©“) since we built the 𝜆.©“/ − 𝜆.©“0 as large as possible, and all equal.

Thus p(1) can be expressed as the initial difference among the smallest pole and the smallest zero (𝜆/“ − 𝜆“0) plus the difference among the other poles and zeroes, that can at most be equal to (𝑛 − 1) ∗ (𝜆/_.©“− 𝜆_.©“0 _{) with 𝜆}

.©“ /

= 1 − 𝜂 and 𝜆.©“0 = “hr_“ − 𝜀, as explained before. In formulas: 𝑝(1) = •𝜆/_“ − 𝜆_“0_{„ + (𝑛 − 1) ∗ (𝜆} .©“ / − 𝜆.©“ 0 _{) = −(}𝑛 − 1 𝑛 + 𝜀) + (𝑛 − 1) ∗ (1 − 𝜂 − 𝑛 − 1 𝑛 − 𝜀) = −𝜀 − (𝑛 − 1)(𝜂 + 𝜀) < 0.

Being 𝑛, 𝜂 𝑎𝑛𝑑 𝜀 all > 0, p(1) will always be negative.

It is straightforward that if the upper limit of the positive values in the equation 𝑝(1) = f 𝜆_./ “ .jr − f 𝜆_.0 “ .jr = (𝜆/_“ − 𝜆0_“_{) + f(𝜆} . / “hr .jr − 𝜆0_.₎

gives a negative p(1), all the other possible pole-zeroes placement in case of •𝜆_“/− 𝜆_“0_{„ >}“hr “ will also cause p(1) to be negative. It is also naïve to understand that if ∑“hr(𝜆_./

.jr − 𝜆0.) gives a negative value (a sum of two negative values gives a negative value).

Linking Proposition 0 to Proposition 2 and 3

(21)

21 a single 0-X block (a pole 𝜆/_. bigger than a zero𝜆.0) big enough and inserted in a combination of X-0 blocks (a combination of pole and zeroes such that ∀𝑖 from 1 to n: 𝜆_./− 𝜆_.0 _{< 0) will turn the} impulse response function (partially) positive (or viceversa, for a X-0 block inserted in a

combination of 0-X blocks).

So, how can we be sure that the impulse response will always be positive? It is straightforward to conclude that by adding 0-X blocks to pole-zeroes combinations of 0-X blocks, might generate an always positive impulse response. Same way backwards, adding X-0 blocks to pole-zeroes

combinations of X-0 blocks would seem logical to produce an impulse response with negative values (at least the first, say p(1)).

Figure 7 - Building ARMA(n,n) through X-0 blocks Figure 8 - Building ARMA(n,) through 0-X blocks

Corollary 0.1: By building up ARMA(n,n) only by adding 0-X to 0-X blocks or X-0 to X-0 blocks (in both cases maintaining the all different poles or all common poles condition), it is possible to obtain an partially positive impulse response (0-X), or to obtain a partially negative impulse response (X-0).

Corollary 0.2: It is immediately clear that by reducing a combination of 0s (zeroes) and Xs (poles) into only 0-X or X-0 couples (that is, 𝒕𝒉𝒂𝒕 ∀𝒊 𝒇𝒓𝒐𝒎 𝟏 𝒕𝒐 𝒏: 𝜆/_. − 𝜆_.0 _{≷ 𝟎, we can prove a partially} positive or partially negative impulse response.

And that’s by trying to extend these considerations and integrating them with Gaalman and Disney (2012, working paper) propositions, here comes Proposition 2 and Proposition 3:

𝐏𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧 𝟐:

Given a ARMA(n,n) model, with n zeroes and n poles, modelling demand as shown by Gaalman (2012, working paper), given no pole-zeroes cancellations (§∀𝑖, 𝑗 𝑓𝑟𝑜𝑚 1 𝑡𝑜 𝑛: 𝜆/_. ≠ 𝜆0»¨) a

sufficient condition to obtain increasing bullwhip as function of the lead time is 𝒕𝒉𝒂𝒕 ∀𝒊 𝒇𝒓𝒐𝒎 𝟏 𝒕𝒐 𝒏: 𝜆/_. − 𝜆_.0 _{> 𝟎.}

(22)

22 𝐏𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧 𝟑:

Given a ARMA(n,n) model, with n zeroes and n poles, modelling demand as shown by Gaalman (2012, working paper), given no pole-zeroes cancellations (§∀𝑖, 𝑗 𝑓𝑟𝑜𝑚 1 𝑡𝑜 𝑛: 𝜆/_. ≠ 𝜆₀_»¨) a

sufficient condition to never obtain increasing bullwhip as function of the lead time is 𝒕𝒉𝒂𝒕 ∀𝒊 𝒇𝒓𝒐𝒎 𝟏 𝒕𝒐 𝒏: 𝜆_𝑖𝜙− 𝜆_𝑖𝜃< 𝟎.

e.g.: X-0 à X-0-X-0 à X-0-X-X-0-0 à X-X-0-X-X-0-0-0 etc.

After having proved proposition 0, and explained the reasoning behind proposition 2 and 3, the following two subsections will show the core analysis of this paper: ARMA(2,2), ARMA(3,3) and some analysis of ARMA(n,n).

ARMA(2,2) and ARMA(3,3). An analytical and simulation perspective

In this subsection we will link proposition 0 to the actual analysis of the Z-Transformed ARMA models. In the previous section we analyzed the link among Proposition 0 and Proposition 2 and 3. Regarding the ARMA(2,2) (with distinct poles) cases enquired by Gaalman and Disney (2009), we will quickly remember that being two poles and two zeroes, we could have 6 possible

combinations, knowing that permutation of 4 elements with 2 repetitions is calculated as follows: 𝑃(𝑛, 𝑖, 𝑗) = 𝑛!

𝑖! ∗ 𝑗!= 4!

2! ∗ 2!= 6.

Of the 6 combinations, 2 behave following proposition 2 (cases 1) and 2), highlighted in green) and 2 behave following proposition 3 (cases 5) and 6), highlighted in red). In particular:

1) 0-X-0-X 2) 0-0-X-X 3) 0-X-X-0 4) X-0-0-X 5) X-0-X-0 6) X-X-0-0

In these cases (ARMA(2,2) with distinct poles) coefficients 𝑟r and 𝑟| are calculated as follows: 𝑟. = ¾•𝑧 − 𝜆./„ ∗ (𝐻(𝑧) − 1)¿ª_‡jˆ ‰ Ž → 𝑓𝑜𝑟 𝑒𝑥𝑎𝑚𝑝𝑙𝑒 𝑓𝑜𝑟 𝑟_r: 𝑟r = ¾•𝑧 − 𝜆r/„ ∗ (𝐻(𝑧) − 1)¿ª_‡jˆ Ž• = •𝑧 − 𝜆r /_{„ ∗}¾•𝑧 − 𝜆r0„•𝑧 − 𝜆0|„ − •𝑧 − 𝜆r/„•𝑧 − 𝜆|/„¿ ¾•𝑧 − 𝜆_r/„•𝑧 − 𝜆/_|„¿ |‡jˆŽ• =¾•𝜆r /_{− 𝜆} r 0_„•𝜆 r /_{− 𝜆} | 0_{„¿ − (𝜆} r /_{− 𝜆} r /_)(𝜆 r /_{− 𝜆} r /₎ •𝜆_r/− 𝜆_|/„

(23)

23 1) 0-0-# sure bullwhip 2) 0-#-0 depends 3) #-0-0 sure no bullwhip

Explicitating coefficient 1 and 2 through the formula: 𝑟.(𝑛) = 1 𝑛 − 𝑖!Á 𝑑“h. 𝑑𝑧“h.•𝑧 − 𝜆/„ “ ∗ (𝐻(𝑧) − 1)ÂÃ ‡jˆŽ And using: 𝑝(𝑘 + 1) = f[𝑟.(𝑛) ∗ “ .jr 𝜆/›h.…r _∗ 𝑘! (𝑘 − 𝑖 + 1)!] Therefore we have (coefficients are generally called with letter “r” or “c”):

𝑟_r(2) = 𝑐_r = •𝜆/_{− 𝜆} r 0_{„ + •𝜆}/_{− 𝜆} | 0_„ 𝑟|(2) = 𝑐| = •𝜆/− 𝜆r0„ ∗ •𝜆/− 𝜆|0„

The discrete time version of the impulse response re-transformed function is:

Which in continuous time gives: 𝑝(𝜏 + 1) = 𝜆/Å∗ 𝑟r(2) + 𝜆/Åhr∗ 𝑟|(2) ∗ 𝜏.

Case 1 - 0-0-#

In this evenience, both r1(2) and r2(2) are positive: hence, p(1) is positive, and will smoothly decrease, ending up at a positive value very close to zero (0+_{). Hence, in this case we have sure} bullwhip. Therefore, it clearly follows Proposition 2.

Case 2 - 0-#-0

In this case we have an always negative r2(2), which makes impulse response always negative for large k (or large t). However r1(2) could either be positive or negative.

Case 3 - #-0-0

The last case sets r1(2) always negative, but r2(2) brings the impulse response positive for large k (or large t). It follows Proposition 3.

(24)

24 Here ARMA(3,3) case will be discussed: a rapid overview of the ARMA(3,3) in case of distinct poles will first be provided, with the help of some simulated plots. An analyitical discussion of

ARMA(3,3) in case of all common poles will follow, and eventually the case of ARMA(3,3) with 2 common poles and a distinct one.

ARMA(3,3) with all distinct poles

In case of all distinct poles we have _}!∗}!Ê! = 20 different dispositions. Of these, 4 should follow proposition 2, and 4 proposition 3.

Of the first 10 disposition that start with a pole (say, the biggest element on the 0-1 line is a pole), only 5 follow proposition 1 (highlighted in green)

1) 0-0-0-X-X-X 2) 0-0-X-0-X-X 3) 0-X-0-0-X-X 4) 0-0-X-X-0-X 5) 0-X-0-X-0-X 6) X-0-0-0-X-X 7) X-0-0-X-0-X 8) 0-X-X-0-0-X 9) X-0-X-0-0-X 10) X-X-0-0-0-X

Of the 10 dispositions that start with a zero (say, the biggest element on the 0-1 line is a zero), only 5 follow proposition 3 (highlighted in red).

11) X-X-X-0-0-0 12) X-X-0-X-0-0 13) X-0-X-X-0-0 14) X-X-0-0-X-0 15) X-0-X-0-X-0 16) 0-X-X-X-0-0 17) 0-X-X-0-X-0 18) X-0-0-X-X-0 19) 0-X-0-X-X-0 20) 0-0-X-X-X-0

Regarding cases 1 to 5, analytical proof of these combinations are hardly possible, since the coefficient calculations often show a negative component whose effect is not clear: from the simulations, however, turns out that all the 5 cases (highlighted in green) always follow proposition 2. The coefficient calculations is made, as always, this way:

𝑟_. = ¾•𝑧 − 𝜆/_. „ ∗ (𝐻(𝑧) − 1)¿ª_‡jˆ ‰ Ž

(25)

25 (𝜆_/r) is always multiplied by a positive coefficient (𝑟r), being the differences among it and the zeroes, and the difference among it and the other poles always positive.

Therefore, simulation comes at a hand: the code A.Y is built implementing the conditions explained in proposition 2, say, having for every i, each pole 𝜆_/. bigger than the correspondent zero 𝜆_0., by these two instructions:

for i = 1:n-1 while P(1,i+1)>=P(1,i) P(1,i+1)=randi([1 99],1); end; end; […] for i = 1:n while Z(1,i)>=P(1,i) Z(1,i)=randi([1 99],1); end; end;

These 2 instructions select only random values which respect proposition 2: having taken 3 poles and 3 zeroes n=3, the first inequality forces poles 2 and 3 to be smaller than pole 1 (P(1,n) with n=3, is the vector of poles). The second inequality, makes the random generation of the zeroes values bound by the limit value of the pole with the same index, say, the first zero has to be smaller than the first pole, and so on, but no limit conditions have been placed regarding each

zeroes and the poles with different index: thus, each of the 5 combinations has been randomly

generated and tested. The only (rare) cases in which any of the zeroes was exactly equal to any of the poles have are not contemplated, thank to the ‘>=’ condition, since they bring the coefficient of the pole 𝜆/_. to 0 (but this has been included in the formulation of Proposition 2). The simulated graphs were checked by the program through the following condition:

for i=1:timespan

if tdata(i)<0 l=1;

end;

(26)

26

Figure 9 - 100 ARMA(3,3) following Prop. 2

Cases 6) to 10) have ambivalent behavior: they may either start with a positive p(1) or not, depending on the pole-zeroes values. Each of them presents a X-0 block, which means that:

− 𝜆_•0_{) , p(1) < 0. Here are provided 20 simulations for each of the} combinations 6) to 10).

(27)

27

Figure 12 - Case 8) Figure 13 - Case 9)

To analyze cases like these, specific conditions had to be implemented (for instance for case 10): P(1,1)=randi([1 99…9],1); Z(1,1)=randi([1 (P(1,1)-1)],1); Z(1,2)=randi([1 (Z(1,1)-1)],1); Z(1,3)=randi([1 Z(1,2)-1],1); P(1,2)=randi([1 (Z(1,3)-1)],1); P(1,3)=randi([1 Z(1,3)-2],1); poles=P/10…0; zeroes=Z/10…0;

(with of course “…” not present in the code, but just to indicate that the level of accuracy can vary).

Regarding cases 11) to 15), being that ∀i from 1 to n: 𝜆_𝑖𝜙− 𝜆𝑖𝜃< 0, it is self evident that they all start from a negative p(1). Therefore, they always follow Proposition 3. In graphs:

Coding the condition ∀ 𝑖 = 1: (𝑛 = 3) ∶ 𝜆_./ −𝜆_𝑖𝜃< 0 only took reverting an inequality in the second while-cycle (the one which determines the position of the zeroes):

for i = 1:n; while Z(1,i)<=P(1,i) Z(1,i)=randi([5 P(1,i)],1); end; end;

Cases 16) to 20) behave mirrorfold to case 6) to 10): in fact, every case from 16) to 20) presents at least a 0-X block, which means: 𝑓𝑜𝑟 𝑖 = 1: (𝑛 = 3) ∃ (𝜆_./,𝜆_.0_) |_𝜆

. /₋_𝜆

.

0 _{> 0.}_{As shown before, p(1) may either be positive or negative.} Simulations for cases 16) to 20) are omitted since they look exactly as the ones of cases 6) to 10). ARMA(3,3) with all common poles

Figure 14 - Case 10)

(28)

28 Let’s now move to the ARMA(3,3) with all common poles. This time, we will face the problem analytically. In fact, the calculation of coefficients changes into a more complex expression: derivation and nonlinearity come into play. In formulas:

𝑟.(𝑛) = 1 𝑛 − 𝑖!Á 𝑑“h. 𝑑𝑧“h.•𝑧 − 𝜆/„ “ ∗ (𝐻(𝑧) − 1)ÂÃ ‡jˆŽ After calculations, we have these explicit expression for 𝑟r(3) , 𝑟|(3) and 𝑟}(3) :

𝑟r(3) = f•𝜆/− 𝜆0.„ “j}

.jr

𝑟_r(3) is logically the sum of poles minus the sum of zeroes, as we can understand from the double derivative (n=3, i=1 à n-i=2) of r_|∗ ∏’j}•𝑧 − 𝜆0_._„

.jr which logically exclude all the polynomial terms

in z with a degree less than 2. In formulas: 𝑑| 𝑑𝑧| 1 2∗ Ì•𝑧 − 𝜆.0„ ’j} .jr Í ‡jˆŽ = 𝑑| 𝑑𝑧|Î 1 2∗ 𝑧}− 𝑧|∗ f •𝜆0.„ ’j} .jr – 𝑧 ∗ … Ð ‡jˆŽ = 3𝜆_/− f •𝜆0_._„ ’j} .jr = f•𝜆/_{− 𝜆} . 0_„ “j} .jr = 𝑟r(3)

This can be extended to every ARMA(n,n) with all common poles, as we will see later on. Below the calculation of the other 2 coefficients (for 𝑟}(3) is pretty straightforward being no derivation at all since n=3 and i=3, n-i=0).

𝑟_|(3) = 𝑟_|(2) + 𝑟_r(2) ∗ (𝜆/_{− 𝜆} } 0_{) = f(Ì•𝜆}/_{− 𝜆} . 0_„ ’j} •jr •©. _. “j} .jr 𝑟_}(3) = Ì•𝜆/_{− 𝜆} . 0_„ ’j} •jr •©.

Therefore in case of 3 common poles, we have the following possible orderings (the 3 common

poles are indicated with #):

1) 0-0-0-# 2) 0-0-#-0 3) 0-#-0-0 4) #-0-0-0

The first case should behave following proposition 2. We will analytically prove that it always does. Therefore the first case confirms increasing bullwhip behavior as a function of lead time.

The last case does never: it will always start from a negative value, invalidating the possibility to have an always positive impulse response.

(29)

29 Beside the pole-zeroes disposition, all the combinations share (obviously) the same time-domain re-transformation: 𝑝(𝑘 + 1) = f[𝑟_.(𝑛) ∗ “ .jr 𝜆/›h.…r _∗ 𝑘! (𝑘 − 𝑖 + 1)!] That after substitutions:

𝑝(𝑘 + 1) = 𝜆/›_{∗ Ñ𝑟}

r(3) + 𝑟|(3) ∗ (𝑘) ∗ 𝜆/hr+ 𝑟}(3) ∗

(𝑘)(𝑘 − 1)

2 ∗ 𝜆/h|Ò

Definition 1.0 still evidently apply, since for k=0 (p(1)) only 𝑟r(3) doesn’t multiply to 0.

Let’s shift from the discrete time version with time variable k to the continuos one, with time variable t. In this case we will call continuous time t as x, for analogy with the general variable utilized for polynomial analysis:

𝑝(𝑥 + 1) = 𝜆/ Ó_{∗ Ñ𝑟} r(3) + 𝑟|(3) ∗ (𝑥) ∗ 𝜆/hr + 𝑟}(3) ∗ (𝑥)(𝑥 − 1) 2 ∗ 𝜆/h|Ò > 0 𝑝(𝑥 + 1) = 𝜆/ Ó_{∗ Ñ𝑟} r(3) + [𝑟|(3)𝜆/hr− 𝑟_}(3) 2 𝜆/h|] ∗ (𝑥) + 𝑟}(3) ∗ 𝑥| 2 ∗ 𝜆/h|Ò > 0 The calculations are hereby developed adopting the following notations:

𝑎𝑥|_{+ 𝑏𝑥 + 𝑐 = 0 => 𝑎 = 𝑟} }(3) ∗ 𝜆/h| 2 , 𝑏 = 𝑟|(3)𝜆/hr− 𝑟_}(3) 2 𝜆/h|, 𝑐 = 𝑟r(3) A second degree complete polynomial in x, to be always positive, requires at least both terms a and c to be positive, being a the dominant factor for x->infinite, and c for t->0. Thus:

𝑐 > 0 → 𝑟_r(3) > 0 => •𝜆/_{− 𝜆} r 0_{„ + •𝜆}/_{− 𝜆} | 0_{„ + •𝜆}/_{− 𝜆} } 0_{„ > 0 → f 𝜆}/ } .jr − f 𝜆_.0 } .jr 𝑎 > 0 → 𝑟}(3) > 0 → Ì •𝜆/− 𝜆.0„ } .jr

(30)

30 𝑥_r,| = −(𝜆/r_{∗ f} 1 𝜆/_{− 𝜆} . 0 “j} .jr − 1 2) ± Û{Á𝜆/r_{∗ ∑} 1 𝜆/_{− 𝜆} . 0 “j} .jr − 12Â | − ∑ ( 2 •𝜆/_{− 𝜆} . 0_„•𝜆/_{− 𝜆} •0„ “j} .jr •©. Ø 2 ∗ 𝑟_}(3) ∗ 𝜆/h|_{∗ 12} Hence we get: 𝑥r,| = −(f 𝜆/ 𝜆/_{− 𝜆} . 0 “j} .jr − 1 2) ± ß{Áf 𝜆/ 𝜆/_{− 𝜆} . 0 “j} .jr − 1 2Â | − f • 2𝜆/| •𝜆/_{− 𝜆} . 0_„•𝜆/_{− 𝜆} •0„ ‘ “j} .jr •©. Ø

So the conditions to have a positive term inside the first parenthesis of the expression above and real value under the square root are:

C1: à∑ _ˆ_Žˆ_hˆŽ

‰ Š “

.jr −r_|á > 0 & if C1 is verified we need to check C2: ∑ â |ˆ

Verifying C1 for case 1 (0-0-0-#) is very straightforward, because 0 < 𝜆/_{< 1 and 0 < 𝜆}

0‰ < 𝜆/, so that 2 ∑ _ˆ_Žˆ_hˆŽ ‰ Š “ .jr > 1 and ∑ ã ˆ ŽØ •ˆŽ_hˆ ‰ Š_„§ˆŽ_hˆ » Š_¨å “j} .jr •©. > 0 , since ˆŽ ˆŽ_hˆ ‰ Š > 1, by construction. And therefore the zeroes could only be found for negative value of t, which are out of our interest. Also r1(3) (=c) is always positive.

Figure 16 - ARMA(3,3) all common poles, case 1)

Going on, we can see immediately that for case 2 (0-0-#-0) the component a (say, r3(3)) is

negative. In fact, r3(3)= ∏ •𝜆}.jr /− 𝜆.0„ is always a product of 2 positive quantities and a negative one (𝜆/_{− 𝜆}

r

(31)

31

Case 4 (#-0-0-0) has both a (r3(3)= ∏ •𝜆}.jr /− 𝜆.0„) and c (𝑟r(3) = ∑}.jr𝜆/− ∑}.jr𝜆.0) always negative.

The most interesting case is therefore Case 3 (0-#-0-0): the fact that a is positive and c might be positive as well, may deceive us: in case of positive c, there might be bullwhip. I originally postulated that there could always be 2 zeroes, since the first simulations provided only results heading that way. However, I mathematically proved that there could be a positive or a negative zero, depending on the values of p(1) (say, 𝑟r(3)).

𝐶1 = ãf 𝜆/ 𝜆/_{− 𝜆} . 0 “j} .jr − 1 2å = 𝜆/ 𝜆/_{− 𝜆} r 0+ 𝜆/ 𝜆/_{− 𝜆} | 0+ 𝜆/ 𝜆/_{− 𝜆} } 0− 1 2 > 0 −→ Case 3): 𝑛𝑜𝑡 𝑎𝑙𝑤𝑎𝑦𝑠 𝐶2 = f • 2𝜆/| •𝜆/_{− 𝜆} . 0_„•𝜆/_{− 𝜆} •0„ ‘ > 0−→ Case 3) ∶ 𝑖𝑓 𝑟_r(3) > 0 “j} .jr •©. E.g.: 𝜆/ _{= 0.48, 𝜆} r 0 _{= 0.55, 𝜆} | 0 _{= 0.51, 𝜆} } 0 _{= 0.03−→ 𝐶1 = −22.295, 𝐶2 = 170.6667}

C2 depends on the value of 𝑟r(3).

Proof: C2 > 0 always for case 3)

f ? â 2𝜆/| •𝜆/_{− 𝜆} 0‰„ §𝜆/− 𝜆0»¨ æ > 0 𝑡ℎ𝑎𝑡 𝑖𝑠 2𝜆 /|_§•𝜆/_{− 𝜆} r 0_{„ + •𝜆}/_{− 𝜆} | 0_{„ + •𝜆}/_{− 𝜆} } 0_„¨ •𝜆/_{− 𝜆} r 0_„•𝜆/_{− 𝜆} | 0_„•𝜆/_{− 𝜆} } 0_„ > 0 𝑡ℎ𝑒𝑟𝑒𝑓𝑜𝑟𝑒 (+) ∗ (𝑟r(3)) (+) ∗ (−) ∗ (−) = 𝑟r(3) > 0 “j} .jr

(32)

32 Hence, in the example, by finding a negative C1, and a positive C2, we calculate the zeroes of the

impulse response, which happen at time 𝑡r = −35.3544 𝑎𝑛𝑑 𝑡| = −9.2266.

Therefore, there may for case 3), there may be increasing bullwhip behavior as a function of lead time in case 𝑟r(3) > 0.

The ARMA(3,3) case: 2 common poles and a single pole.

The case in which an ARMA(3,3) model has a double pole and a single pole is controversial in many ways. The “all common poles”-case always had no differentiating variable at the denominator: this feature makes inference over the bullwhip a bit more difficult, because it adds a negative

component to the residuals calculation. Furthermore, the calculations become tricky but still proposition 2 or 3 are applicable because p(1) is equal to the sum of poles minus the sum of zeroes.

Proposition 2 and 3 seems to hold for ARMA(3,3) with 2 common poles. The thousands simulation launched never reported a negative value of the impulse response for the cases that (should) theoretically follow Proposition 2. Proposition 3 has never been contradicted either.

In the same way, other inferences can be done through simulation for ARMA(n,n) with n>3 with 1<j<n common poles. Let’s proceed first expliciting p(t+1) in t=0, say p(1): the residual calculation is performed below.

The notation to indicate the common poles is the following: r1,2(3) for the residual of the first of the two common poles (the derivated one), r2,2(3) for the second common pole residual (the non derivated one) and r3 for the non common, third, pole residual. The two common poles will be indicated as 𝜆/_{and the third pole as 𝜆}

(33)

33 𝑟r,|(3) = 𝑑 𝑑𝑧Î•𝑧 − 𝜆/„ |•𝑧 − 𝜆r0„•𝑧 − 𝜆0|„•𝑧 − 𝜆0}„ (𝑧 − 𝜆/₎|_{•𝑧 − 𝜆} } /_„ Ð |‡jˆŽ = 𝑑 𝑑𝜆/Î•𝑧 − 𝜆/„ |•𝜆/− 𝜆_r0„•𝜆/− 𝜆_|0„•𝜆/− 𝜆_}0„ •𝑧 − 𝜆_/„|•𝜆/_{− 𝜆} } /_„ Ð |‡jˆŽ = 𝑑 𝑑𝜆/Î 𝑟}(3) •𝜆/_{− 𝜆} } /_„Ð = Î𝑟|(3) ∗ •𝜆 /_{− 𝜆} /í„ − 𝑟}(3) •𝜆/_{− 𝜆} /í„ | Ð 𝑟_|,|(3) = Á•𝑧 − 𝜆/_„|•𝑧 − 𝜆r0„•𝑧 − 𝜆0|„•𝑧 − 𝜆0}„ (𝑧 − 𝜆/₎|_{•𝑧 − 𝜆} /í„ Â |_‡jˆ_Ž = Î•𝑧 − 𝜆/„| •𝜆/_{− 𝜆} r 0_„•𝜆/_{− 𝜆} | 0_„•𝜆/_{− 𝜆} } 0_„ (𝑧 − 𝜆/₎|_{•𝑧 − 𝜆} } /_„ Ð |‡jˆŽ = Î 𝑟_}(3) •𝜆/_{− 𝜆} } /_„Ð 𝑟_} = Î•𝑧 − 𝜆/_}„|•𝑧 − 𝜆r0„•𝑧 − 𝜆|0„•𝑧 − 𝜆}0„ (𝑧 − 𝜆/₎|_{•𝑧 − 𝜆} } /_„ Ð |‡jˆŽí = Î•𝑧 − 𝜆/_}„•𝜆} /_{− 𝜆} r 0_„•𝜆 } /_{− 𝜆} | 0_„•𝜆 } /_{− 𝜆} } 0_„ (𝑧 − 𝜆/₎|_{•𝑧 − 𝜆} } /_„ Ð |‡jˆ_Ží = Î ∏} •𝜆_/_í− 𝜆_î_‰„ .jr •𝜆/_} − 𝜆/_„ Ð For t=0, p(t+1)=p(1) and therefore:

𝑝(𝑡 + 1) = ï𝑟_r,|(3)𝜆/b_{+ 𝑟}

|,|(3) ∗ 𝑡 ∗ 𝜆/bhr + 𝑟}∗ 𝜆}/bð = 𝑝(1)|bjk = = {𝑟_r,|(3) + 𝑟_}}

A neater expression, which highlights the validity of my statement on p(1) position also for this apparently complex expression, can be obtained by calling 𝜆_}/− 𝜆/_{as 𝜀. In fact:}

(34)

34 𝑝(1) = 𝑟_r(3) + 𝜀 = 𝜆/_{− 𝜆} r 0 _{+ 𝜆}/_{− 𝜆} | 0 _{+ 𝜆}/_{− 𝜆} } 0 _{+ 𝜀 = 𝜆}/_− 𝜆 r 0 _{+ 𝜆}/_− 𝜆 | 0 _{+ (𝜆}/_{+ 𝜀) − 𝜆} } 0 = 𝜆/_− 𝜆 r 0 _{+ 𝜆}/_− 𝜆 | 0_{+ 𝜆} /}− 𝜆}0 = f(𝑝𝑜𝑙𝑒𝑠) − f(𝑧𝑒𝑟𝑜𝑠) } .jr } .jr

Let’s distinguish 20 orderings : !.

With positive 𝜀, that is 𝜆/_} > 𝜆/_:

1) à 0-0-0-#-X 5)à #-X-0-0-0 2)à 0-0-#-0-X 6)à #-0-X-0-0 3)à 0-#-0-0-X 7)à #-0-0-X-0 4)à #-0-0-0-X 8)à 0-#-0-X-0 9)à 0-#-X-0-0 10)à 0-0-#-X-0 With negative 𝜀, that is 𝜆/_} < 𝜆/_:

11)à 0-0-0-X-# 15)à X-#-0-0-0 12)à 0-0-X-0-# 16)à X-0-#-0-0 13)à 0-X-0-0-# 17)à X-0-0-#-0 14)à X-0-0-0-# 18)à 0-X-0-#-0 19)à 0-X-#-0-0 20) à 0-0-X-#-0

Proposition 2 apparently always holds: in fact, following proposition 1 and 2, the impulse response of the green-highlighted (1,2,11,12,13) cases should always be positive: a 10’000 simulated data say so. Figure X shows 100 graphs of combinations (1,2,11,12,13), which follow proposition 2. The plots are build leaving no conditions on pole-zero orderings beside Proposition 2 ones: ∀ 𝑖 = 1: (𝑛 = 3) 𝜆_𝜙 𝑖 − 𝜆𝜃𝑖 > 0. In code: P(1,1)= randi([5 999],1); P(1,2)=P(1,1); P(1,3)=randi([5 999],1); for i = 1:n; while Z(1,i)>=P(1,i) Z(1,i)=randi([1 P(1,i)],1); end; end; poles=P/1000; zeroes=Z/1000;

Despite the hard trying, none of the combinations following Proposition 2 (1,2,11,12,13) ever hit a value below zero. Therefore, it is straightforward to conclude that proposition 2 holds for every ARMA(3,3) model, provided no cancellations occur §∀𝑖, 𝑗 𝑓𝑟𝑜𝑚 1 𝑡𝑜 𝑛: 𝜆_𝑖𝜙≠ 𝜆_0»¨.

(35)

35

Figure 19 - ARMA(3,3) with 2 common poles following Prop.2

It is superfluous to analyze the red highlighted cases (, since they all start from a negative p(1). However, I provided 20 simulated graphs built by implementing the condition of Proposition 3 (∀ 𝑖 = 1: (𝑛 = 3) 𝜆_./− 𝜆_.0 _{< 0)}_{and of course} §∀𝑖, 𝑗 𝑓𝑟𝑜𝑚 1 𝑡𝑜 𝑛: 𝜆_𝑖𝜙≠ 𝜆𝑖𝜃¨. In code: P(1,1)= randi([5 999],1); P(1,2)=P(1,1); P(1,3)=randi([5 999],1); for i = 1:n; while Z(1,i)>=P(1,i) Z(1,i)=randi([1 P(1,i)],1); end; end; poles=P/1000; zeroes=Z/1000;

The cases who do not follow Proposition 2 or 3 (3,4,8,9,10,14,17,18,19,20) do behave exactly as