The Tremor Stochastic-Local-Volatility Model

合集下载

volatility surface

1
risk when these options are priced. Traders also use the volatility surface in an ad hoc way for hedging. They attempt to hedge against potential changes in the volatility surface as well as against changes in the asset price. As described in Derman (1999) one popular approach to hedging against asset price movements is the “volatility-by-strike” or “sticky strike” rule. This assumes that the implied volatility for an option with a given strike price and maturity will be unaﬀected by changes in the underlying asset price. Another popular approach is the “volatility-bymoneyness” or “sticky delta” rule. This assumes that the volatility for a particular maturity depends only on the moneyness (that is, the ratio of the price of an asset to the strike price). The ﬁrst attempts to model the volatility surface were by Rubinstein (1994), Derman and Kani (1994), and Dupire (1994). These authors show how a one-factor model for an asset price, known as the implied volatility function (IVF) model, can be developed so that it is exactly consistent with the current volatility surface. Unfortunately, the evolution of the volatility surface under the IVF model can be unrealistic. The volatility surface given by the model at a future time is liable to be quite diﬀerent from the initial volatility surface. For example, in the case of a foreign currency the initial U-shaped relationship between implied volatility and strike price is liable to evolve to one where the volatility is a monotonic increasing or decreasing function of strike price. Dumas, Fleming, and Whaley (1997) have shown that the IVF model does not capture the dynamics of market prices well. Hull and Suo (2002) have shown that it can be dangerous to use the model for the relative pricing of barrier options and plain vanilla options. In the ﬁrst part of this paper we develop a general diﬀusion model for the evolution of a volatility surface and derive a restriction on the speciﬁcation of the model necessary for it to be a no-arbitrage model. Other researchers that have independently followed a similar approach are Ledoit and Santa Clara (1998), Sch¨ onbucher (1999), and Brace et al (2001). In addition, Britten–Jones and Neuberger (2000) produce some interesting results characterizing the set of all continuous price processes that are consistent with a given set of option prices. Our work is diﬀerent from that of other researchers in that we a) investigate the implications of the no-arbitrage condition for the shapes of the volatility surfaces likely to be observed in diﬀerent situations and b) examine whether the various rules of thumb that have been put forward by traders are consistent with the no-arbitrage condition. We also extend the work of Derman (1999) to examine whether the rules of thumb are supported 2

Collective Transport From Superconductors to Earthquakes

a rXiv:c ond-ma t/9711179v118Nov1997Collective Transport in Random Media:From Superconductors to Earthquakes Lectures at Summer School on “Fundamental Problems in Statistical Mechanics IX”August 20-23,1997Daniel S.Fisher Physics Department,Harvard University Cambridge,MA 02138ﬁsher@ Abstract In these lectures,a variety of non-equilibrium transport phenom-ena are introduced that all involve,in some way,elastic manifolds be-ing driven through random media.A simple class of models is studied focussing on the behavior near to the critical “depinning”force above which persistent motion occurs in these systems.A simple mean ﬁeld theory and a “toy”model of “avalanche”processes are analyzed and used to motivate the general scaling picture found in recent renormal-ization group studies.The general ideas and results are then applied to various systems:sliding charge density waves,critical current be-havior of vortices in superconductors,dynamics of cracks,and simple models of a geological fault.The roles of thermal ﬂuctuations,defects,inertia,and elastic wave propagation are all discussed brieﬂy.I.IntroductionMany phenomena in nature involve transport of material or some other quan-tity from one region of space to another.In some cases transport occurs in systems that are close to equilibrium with the transport representing only a small perturbation such as ﬂow or electrical current in a metal,while in other cases it involves systems that are far from equilibrium such as a land-slide down a mountain,or a drop of water sliding down an irregular surface.Sometimes,particles or other constituents move relatively independently of1each other like the electrons in a metal,while in other situations the inter-actions play an important role,as in the landslide and the water droplet.If the interactions are strong enough,all the particles(or other con-stituents)move together and the macroscopic dynamics involves only a small number of degrees of freedom.This is the case for a small water drop on,e.g. wax paper,which slides around while retaining its shape.But if the interac-tions are not so strong relative to the other forces acting on the constituents, then the transport involves in an essential way many interacting degrees of freedom.This is the case for a larger water drop on an irregular surface for which the contact line between the droplet and the surface continually de-forms and adjusts its shape in response to the competition between the surface tension of the water and the interactions with the substrate[36][16].Such a moving drop and the landslide are examples of non-equilibrium collective transport phenomena,which will be the general subject of these lectures.This is,of course,an impossibly broad subject!We must thus narrow the scope drastically.Although the range of systems discussed here will, nevertheless,be reasonably broad,we will primarily focus on systems in which the interactions are strong enough so that the transported object(or at least some part of it)is elastic.We will use this in a general and somewhat loose sense that the transported object has enough integrity that if one part of it moves a long distance then so,eventually,must the other parts as well. Thus theﬂuid drop is elastic if it does not break up—i.e.its perimeter retains its integrity—while a landslide is not elastic as some rocks will fall much further than others and the relative positions of the rocks will be completely jumbled by the landslide.We will be interested in systems in which the medium in which the transport occurs has static random heterogeneities(“quenched randomness”) which exert forces on the transported object that depend on where it is in space.Examples we will discuss are:interfaces between two phases in random media[28][2],such as between twoﬂuids in a porous medium[37],or do-main walls in a random ferromagnetic alloy;lattices of vortices in dirty type II superconductors[34];charge density waves which are spatially periodic modulations of the electron density that occur in certain solids[35][4][3];and the motion of geological faults[43],[5].In addition to the contact line of theﬂuid drop already mentioned[36],[16],another well known—but poorly understood—example that we will,however,not discuss is solid-on-solid fric-2tion.In all of these systems,a driving force,call it F,can be applied which acts to try to make the object move,but this will be resisted by the random “pinning”forces exerted by the medium or substrate.The primary questions of interest will involve the response of the system to such an applied driving force[47].If F is small,then one might guess that it will not be suﬃcient to overcome the resistance of the pinning forces;sections of the object would just move a bit and it would deform in response to F,but would afterwords be at rest.[Note that in most of what follows,we will ignoreﬂuctuations so that the motion is deterministic and the objects can be said to be stationary.] If the force is increased,some segment might go unstable and move only to be stopped by stronger pinning regions or neighboring segments.But for large enough F,it should be possible to overcome the pinning forces—unless they are so strong that the object is broken up,an issue we will return to at the end—and the object will move,perhaps attaining some steady state velocity v.Basic questions one might ask are:is there a unique,history independent force,F c separating the static from the moving regimes?How does v depend on F(and possibly on history)?Are there some kinds of non-equilibrium critical phenomena when v is small?How does the system respond to an additional time or space dependent applied force?These are all macroscopic properties of the system.But we will also be interested in some microscopic properties:how can one characterize(statistically)the deformations of the object when it is station-ary[38]?The dynamic deformations and local velocities when it is moving? The response to a small local perturbation?etc.Motivated by possible analogies with equilibrium phase transitions[33], we can ask if there are scaling laws that might obtain near a critical force which relate,for example,the characteristic length scale L,for some process, to its characteristic time scale,τ,via a power law relation of the form:τ∼L z(1) Trying to answer some of these questions—and to pose other more pointed questions—is the main aim of these lectures.In the next few sections a particular system and its natural(theoretical)generalizations will be studied and tools and ideas developed.In the last section,these are tentatively3applied to various physical systems and some of the complicating features left out of the simple model systems are discussed.This leads naturally to many open questions.II.Interfaces and ModelsIn order to develop some of the general ideas—both conceptual and compu-tational—we will focus initially on an interface between two phases that is driven by an applied force through an inhomogenous medium[28],[2].The essential ingredients of a model of this system are:the forces of sections of the interface on nearby sections,i.e.the elasticity of the interface caused by its interfacial tension;the preference of the interface for some regions of the system over others due to the random heterogeneities;and some dynamical law which governs the time evolution of the local interface position.We will initially make several simplifying approximations,which we will come back and examine later.First,we assume that the interface is not too distorted away from aﬂat surface normal to the direction(z)of the driving force so that its conﬁguration can be represented by its displacementﬁeld u( r)away from aﬂat reference surface.The coordinates R=(x,y,z)of points on the surface are then(x,y)= r and z=u( r).(2) Second,we will assume that the dynamics are purely dissipative i.e.that iner-tia is negligible—a good approximation in many physical situations.Keeping only the lowest order terms in deviations fromﬂat,we then have∂u( r,t)ηFigure1:Schematic of a one-dimensional interface in a two dimensional disordered system illustrating the forces acting on the interface.the“stress”on the interface from its elasticity which is“transmitted”by the kernel J( r,t).Short range elasticity of the interface corresponds toJαδ(t)∇2δ( r).(5) A schematic of such an interface and the forces acting on it is shown in Fig1.Keeping in mind some of the other problems of interest[36],[6,30]in addition to interfaces,we will abstract to a more general problem of a d-dimensional elastic“manifold”—d=2for the interface—with more general interactions,which can be long-range,embodied in J( r,t).In addition tothe form of J( r,t),the system will be characterized by the statistics of the pinning forces which impede interface motion near points where the interface has lower(free)energy;f p( R)will generally have only short-range correla-tions in space,i.e.,in both u−u′and r− r′.Even with these simplifying assumptions,the model Eqs(3,4)is impossible to analyze fully due to the non-linearities implicit in the u dependence of f p( r,u).Nevertheless a lot ofthe qualitative behavior can be guessed.If the driving force is suﬃciently small,then it will be insuﬃcient toovercome the pinning forces.But if F is increased slowly,it may overcome the pinning of some small segment of the interface which can then jump forwards only to be stopped by stronger pinning forces or by the elastic forces from neighboring still-pinned parts of the interface.But if the drive5is larger,the neighboring regions may themselves not be strongly enough pinned to resist the increase in stress from the jumping section and may themselves jump forward leading to an“avalanche”of some larger region of the interface;this process might or might not eventually stop[38].If the force is large enough—and certainly if it exceeds the maximum f p—then it is not possible for the interface to be pinned and the interface will move forward with some average velocity∂t [u a( r1,t)−u b( r1,t)]|t=t1=σ[ r1,{u a}]−σ[ r1,{u b}]= d r′ t1dt′J( r1− r′,t1−t′)[u a( r′,t′)−u b( r′,t′)]>0(6) since the pinning force at r1is the same in both conﬁgurations and therefore cancels out.By assumption the last expression in Eq(6)is positive as long as J is non-negative so that for t>t1,u a is again ahead of u b violating the assumption.The condition thatJ( r,t)≥0(7) for all r,t plays an important role in the theoretical analysis and frequently also in the physics of these types of systems.We will refer to models with this convexity property as monotonic;they have the property that if the displacements and the driving force F(t)increase monotonically with time, then so will the total“pulling force”—see later—on any segment.Except in theﬁnal section we will focus solely on monotonic models.We have shown that in monotonic models one conﬁguration that is ini-tially behind cannot“pass”another that is ahead of it[26];therefore station-6ary and continually moving solutions cannot coexist at the same F;therefore F c is unique.This is a big simpliﬁcation and one that will not occur gen-erally,in particular not in some of the systems that we discuss in the last section.For forces well above F c,one can use perturbative methods to study the eﬀects of the random pinning and compute,for example:the mean velocity,to u( r)there can be multiple values of u which satisfy Eq(10);as we shall see these play an important role in the physics.At this point,it is helpful to be more concrete.Let us consider a simple model of the pinning consisting of pinning sites u pα( r)distributed forﬁxed r with random spacings between the u pα( r),Υα( r)≡u pα+1( r)−u pα( r)(11) drawn,for each r,independently from a distributionΠ(Υ)dΥ.The pinning force f p[( r,u( r)]=0except if u( r)is equal to one of the pinning positions, while for u( r)=u pα( r),f p can take any value between zero and a yield strength,f y,which is the same for each pin.A typical realization of the pinning force f p(u)on some segment of the interface is plotted in Fig.2a. Note that for aﬁxedφ,there are several possible values of u given by the intersection of the lineφ−˜Ju with f p(u).Ifφis increased,then the particular(history dependent)force-balanced position u(φ)that the interface point is following adiabatically can become unstable—for example,the conﬁguration denoted by the circle in Fig.2a—and u must jump to a new position.During the jumpη∂u∂t =0in the adiabatic parts for the pinning model illustratedin Fig2),is to requireu[φ( r,t)]=u ad[φ( r,t−t d)](12) with someﬁxed(microscopic)delay time t d.This is illustrated in Fig2a. Note that,formally,this can be accomplished,by takingη→0and J( r,t)= J( r)δ(t−t d).III.Inﬁnite-range model:meanﬁeld theoryThe above discussion in terms of the local pulling forceφ( r,t)suggests that we could try to analyze the system crudely by assuming that the spatial and temporalﬂuctuations inφ( r,t)are small so thatφcan be replaced by8fi-fuiFigure2:a)Simple model of the forces on one segment of an interface. The segment can be pinned at the positions,u p iα,of the vertical lines at which the pinning force can take any value up to the yield strength f y.The intersections of the“comb”representing the pinning force f p(u i)and the diagonal lineϕ−˜Ju i withϕthe total pulling force from the applied force and other segments of the interface,are the possible stationary positions of u i indicated by the dots.The one of these with the smallest u i,u m i,plays a special role as discussed in the text.The amount∆ϕthatφneeds to increase by to depin the segment from this pinning position is w i˜J.b)Dynamics of the same segment of the interface as the pulling force is increased.The actual u i(t)(dotted),the adiabatic approximation to this(solid),and the time delayed approximation(dashed)that is used in the analysis in the text are all shown.9some sort of time dependent average(13)Ni.e.,inﬁnite range forces.(Note that Eq(13)includes a self-coupling piece but its eﬀects are negligible in the desired N→∞limit.)Much can be done for general non-negative J(t)and more complicated forms of f p(u)using the actual dynamical evolution Eq(3)[3],but to keep things simple we will use the time-delayed adiabatic approximation discussed above withJ(t)=˜Jδ(t−t d);(14) and the form of f i p of Fig2a with independent randomness for each i.For simplicity,we will focus on the strong pinning limit which corresponds tof y>˜JΥmax.(15) It is left to the reader to show that including some of the more“realistic”features within the inﬁnite range model does not change the qualitative or other universal aspects of the results.10Our task is now simple,at least in principle:we assume some meanﬁeldN i u i(t−t d),(16) compute the evolution of each u i(t)(fromφi=φ(t)until the computed<u(t)>≡1φ(t+t d)−F /˜J for all t.Weﬁrst try the simplest possibility:a constantφ−F)/˜J.(18) But we must be careful:If we start choosing too many large u i’s,we mayﬁnd thatwill become too large.We can thus ask:what are the minimum and maximum possiblefor a givenφ) corresponds to theﬁrst pinning position—i.e.one of the{u p iα}—to the right of the intersection of the line f=φ−˜Ju with the line f=f y that passes through the tips of the“comb”–representing the yield strength–in Fig.2a1. Since the peaks are randomly positioned,min=(19) with w i the distance to the next pin which has the distribution2:Prob(w)= ∞w11The strong pinning condition f y>˜JΥmax ensures that u mi is at a pinning position.The general case can be worked out similarly2We use notations like“Prob(w)”to mean the probability that the continuous variable w is in the range w to w+dw,divided by dw;i.e.Prob(w)is the probability density(usually called by physicists“distribution”)of w.One must remember,however,that if variables are changed e.g.from w to w′,then there is a Jacobian needed:Prob(w′)= dwHere the quantity in parentheses is the probability distribution that a random point is in an interval of widthΥbetween pins;this includes the factor of Υ/Υ≡ ∞0ΥΠ(Υ)dΥ(21) because of the presence of more points in wider intervals.Integration of Eq(20)by parts yields<w>=Υ)so thatmin=(Υ2Υ=F˜J+2Υ2Υ(23)a non-trivial result for the critical force above which no static solutions are possible.Note that as the interaction strength,˜J,is increased,the critical force decreases.Physically,this is a consequence of the elasticity causing the system to average over the randomness more eﬀectively:pulling a stiﬀobject over a rough surface is easier than pulling aﬂexible one.For F<F c the number of stable solutions,N s will be exponentially large with an“entropy”per segment ln N sφincreasing,we can understand the behavior from Fig2.Asφ).Furthermore,after a jump u will again be on the new smallest u p iαfor the increasedφ)to the next. But now the time delays must play a role.If we assume a solution which progresses uniformly on average,=vt,then[Note that to ensure that u does not stop between pins,we again need the strong pinning condition˜JΥmax<f y]With all u i=u m i(F−F cv=v∼(F−F c)β(26) with the critical exponentβ=βMF=1(27) in this inﬁnite range mean-ﬁeld model[41].Note that a comparison of Eq(25)for F>>F c and the original dynamic equation(3),suggests that a natural choice is t d=ηv=Fv not strictly linear for F>F c.Nevertheless near F c,the mean velocity will still be characterized by the exponentβ=1.A typical mean-ﬁeldVcFigure3:Velocity versus driving force in a typical meanﬁeld model is indi-cated by the solid line.Note the linear dependence of v on F just above F c. The dashed line is the behavior in the absence of pinning.Figure4:Schematic of hysteresis loops that occur as the force is increased from zero to the critical force,decreased to the critical force in the opposite direction,and then cycled between these values.The direction of change of F is indicated by the arrows with the“1”denoting theﬁrst increase.14IV.Avalanche statistics and dynamicsIn the meanﬁeld model introduced in the previous section,the statistics and other properties of the avalanches of jumps that occur as the driving force is increased slowly can be worked out in substantial detail[39][44].We will carry out the analysis using methods which can be generalized to provide useful information about the behavior with more realistic interactions.Let us consider what happens for F<F c when F is increased by a very small amount.If the increase is suﬃciently small,then no segments will jump.But a slightly bigger increase—typically of order1φwill have only advanced by an amount of order 1φis increased by a small amount∆ϕi.¿From Fig.2a,we see that∆ϕi=˜Jw i.For large N,all but very special ways of preparing the conditions before the avalanche starts will yield a distribution of these small∆ϕi which are independent and randomly distributed with(initial-condition dependent)densityρ≡ρ(∆ϕi=0);(28)ρthus measures a local susceptibiltiy to jumping.We can now immediately conclude something about the mean number of jumps<n t>at a time t after the initial jump.Since the n t−1jumps at time t−1will cause an increase inΥsince the distribution ofΥ’s for the almost unstable segments and hence <Υ>could depend on the initial conditions).This will cause,on average,ρ<Υ>n t−1jumps at time t,i.e.<n t>=ρ<Υ>˜J<n t−1>.(29) The crucial parameter is thusρ<Υ>˜J;if this is greater than one the avalanche will runaway.If the system is below F c as we have assumed,it will eventually be stopped only when aﬁnite fraction of the segments have15jumped and the system has found a stable—and more typical—conﬁguration. But ifρ<Υ>˜J<1,then the expected total size of an avalanches≡ i∆u i(30)is simply<Υ><s>=contains the information of interest.Note that vertical bars as in Eq(37) denote“given”;i.e.conditional probability.It is useful to deﬁne a generating function of the distribution including all times up to TΓT{µt}≡<exp(i T t=0µt m t)>P;(38)(usually we will drop the P).Note thatΓT is simply the Fourier transform of P restricted to times≤T.We can derive a recursion relation forΓT in terms ofΓT−1.For a given m T−1,the number of jumps triggered at time T will be Poisson distributed with meanρm T−1i.e.Prob(n T|m T−1)=e−ρm T−1(since n0=1).All the information has thus gone intoλ0.As long as the system is stable,i.e.ρ≤1,we can simply take T→∞to recover the full information.[Ifρ>1,then there is a non-zero(and computable)possibility that s=∞,and more care is needed.]To get the probability distribution of s,we simply set allµt=µand thenProb(s)= µe−iµs e iλ∗(µ)(45) with µ≡1−ǫ2+2ibµ2iǫ2/b for large s and we thushave,after replacing the dummy variableµbyµ−12iµ/b e−iµs e−sǫ2/(2b)(50)By“power counting”,we see that the branch cut must yield a12depen-dence;hence for large s,Prob(s)∼e−sǫ2/(2b)2(51)18[39].Note that for smallǫ,the mean<s>is dominated by large s avalanches which are rare sinceProb(s∼1ǫ2)∼ǫ.(52)We see from Eq(51)that these yield<s>∼1s!,(53) with,of course,P(s=0)=0.From the limiting large s forms!≈s s e−s√ρ(ρe−ρ)s factor in Eq(53).It is nice that theexact result can be found in this case,but in general,asymptotic methods like those we have used above give more understanding and are more widely applicable.Nevertheless,to convince skeptical colleagues,a few exact results are useful!In addition to the distribution of avalanche sizes,we are also interested in their temporal evolution.For example,one might ask what is<m t|s>, i.e.what is the time development of an average event of size s?This can be computed using the generating function.If we chooseµτ=µforτ=t(55)andµt =µ+νt ,(56)then 1∂νt νt =0=<n t e iµs >=∂λ0∂νt = ∂λt ∂λt µt −1×...× ∂λ0∂νtv t =0= ρe iλ∗(µ) t (60)After shifting µas in Eq(50)we see that,for ǫsmall and s and t large,<m t |s >∼12iµ.(61)This will be dominated by µ∼1s.(62)Note that this is much less then the maximum possible duration τmax =s −1.The integral in Eq (61)can be done exactly (by writing µ=−i x22s (63)for large s independent of ρ.Again,the behavior for large s and 1<<t <<s is generic up to a coeﬃcient b that should appear as in Eq.(51).For the20particular constantΥcase,the exact result can be computed from Eq(58) yielding<m t|s>=t+1√depend on the past history.But on a generic approach to F c from below(e.g. after“training”the system by a slow increase to F c from F=−F c),ρwill approach unity at F c and the cutoﬀ˜s∼1.(66)1−ρΥJ(q,ω)The critical point is thus still given by(ρΥ˜J)crit=1(67)22with˜J=J(q=0,ω=0)(68) and the mean total size is<s>=ΥP rob(s) µq ωe−iµs e iλ∗(µ)e i q· r e−iωt2πs3e−isµ12iµ+K( q,ω)= dλe−1s3λ2+K2( q,ω) (72)For an interface with dissipative dynamics and local elasticity,in the absence of pinning or driving forces we have,after rescaling lengths and times,simply∂uBut to understand the general behavior,and to apply the results to other physical systems,we would like to include the possibility of long range elas-ticity,i.e.,dtJ( r,t)∼1(78)[K s(q)]2which is of order one independent of s ifd>d c(α)=2α.(79) We thus see the appearance of a special critical dimension above which no segment will jump more than a few times even in an arbitrarily large avalanche.Indeed,we will see that above the critical dimension driven in-terfaces will have only bounded small scale roughness.For d<d c,the integral in Eq(78)is inﬁnite so that small q(i.e.small K)dominates and more care is needed.The cutoﬀof q of Eq(72)whenand hence q∼s−1K∼λyields,with a typicalλ∼1s2α.(80) Note the appearance of a non-trivial exponent relating∆u and s.It depends, as is usually the case for critical phenomena,on the spatial dimension.As24mentioned in the Introduction,Eq(80)is just the kind of scaling law we expect near critical points.Equation(80)relates characteristic scales of displacement to the characteristic scales of avalanche size.We can also say something about the spatial extent and shape of large avalanches,by computing<∆u(r)|s>.For d>d c(α),the integral in Eq(72)will be cutoﬀfor q>r by the e i q· r oscillations and hence dominated by q∼1r d−2α(81)for1<<r<<s1r )∼λ∼1s.For d<d c,on the other hand,as long as r<<s1λ2+K2sfor typical λwill be dominated by q∼s−12σfor r<<s12α(83) is thus some measure of the diameter of an avalanche.Let us now try to interpret these results[Note that the skeptic could compute e.g.,<[∆u(r)]2|s>etc.to provide further support for the picture below].For d>d c,the fact that<∆u(r)>is much less than unity for r>>1strongly suggests that most segments will not jump even if they are within r<<L of the avalanche center,rather only a fraction∼1/r d−2αof them will jump,and these typically only once or a few times.The number of sites that have jumped at all within a distance R<L of the origin is of order R2α<<R d so that the avalanche is fractal.The total number of sites that jump,its“area”A is thus,by taking R∼L,A∼L d f∼s<<L d(84) with the fractal dimensiond f=2αfor d>d c.(85)In lower dimensions,the picture is quite diﬀerent.The approximate in-dependence of<∆u(r)|s>of r for r<<L suggests that each site in this25region jumps a comparable number of times∼s1−d/(2α)(withﬂuctuations around this of the same order)and hence the avalanche is not fractal but has areaA∼L d∼s d(89)Lκwithκ=α.(90) The duration of an avalanche with dissipative dynamics—corresponding toK(q,ω)≈−iω+|q|σ(91)—is given simply by scaling,i.e.,τ∼L z∼s1s is(not surprisingly)the same as in the inﬁnite-range mean-ﬁeld model.We have found that in our toy model,many of the properties of large avalanches near to the critical point(actually any large avalanche although they are rare away from criticality),obey scaling laws which relate various characteristic physical properties to each other by power law relationships. For example,for d<d c,an avalanche of diameter L has typical size s∼L2α/d,26displacement∆u∼Lζand duration L z.This type of scaling behavior is one of the key aspects of critical phenomena in both equilibrium and non-equilibrium systems.But there is more:if we scale all lengths by a correlation lengthξ∼ǫ−1/α(94) which is the diameter above which avalanches become exponentially rare,andcorrespondingly displacements byξζ,durations byξz,etc.,thenfunctionssuch as those that occur in the distribution of avalanche sizes Eq(51),or the average growth of the displacements during an avalanche,<∂u( r,t)∂t|s>≈C uξ,tC uξd+ζ (95)with s= d r∆u( r)the total size,C u and C t non-universal(dimensionfull) coeﬃcients which set the scales of the displacements and times;these depend on the random pinning,η,D;etc.The universal scaling function isY( R,T,m)= Q Ω dΛe i Q· R−iΩT e−1m3Λ2+(−iΩ+|Q|α)2 (96) which depends only on the dimension,the range of interactions,and the type of dynamics(i.e.dissipative),as is manifested in the low frequency form of the stress transfer function K( q,ω).As we shall see in the next section,a similar scaling structure is expected to exist in more realistic models.Let us now try applying the toy model results to the interface problem with d=2and short range elasticity,so thatα=2.This dimension is less thand short−rangec=d c(α=2)=4,(97) so we haveζ=2,i.e.∆u(L)>>L,for large avalanches.But this is clearly unphysical:our original model for the interface assumed that it was close to ﬂat so that,at least on large scales,we need small angles of the interface i.e.∇u<<1.Thus the result Eq(87)violates the assumptions of our original model in this case.What has gone wrong?Is the original model bad or have we made some grievous errors in trying to analyze it?The answer is the latter and under-standing why gives some clues as to how to do better.27。

The Detection and Characterization of cm Radio Continuum Emission from the Low-mass Protost

a r X i v :0705.1747v 1 [a s t r o -p h ] 12 M a y 2007Draft version February 1,2008Preprint typeset using L A T E X style emulateapj v.03/07/07THE DETECTION AND CHARACTERIZATION OF CM RADIO CONTINUUM EMISSION FROM THELOW-MASS PROTOSTAR L1014-IRSYancy L.Shirley 1,Mark J.Claussen 2,Tyler M.Bourke 3,Chadwick H.Young 4,Geoffrey A.Blake 5Draft version February 1,2008ABSTRACTObservations by the Cores to Disk Legacy Team with the Spitzer Space Telescope have identiﬁed a low luminosity,mid-infrared source within the dense core,Lynds 1014,which was previously thought to harbor no internal source.Followup near-infrared and submillimeter interferometric observations have conﬁrmed the protostellar nature of this source by detecting scattered light from an outﬂow cavity and a weak molecular outﬂow.In this paper,we report the detection of cm continuum emission with the VLA.The emission is characterized by a quiescent,unresolved 90µJy 6cm source within 0.′′2of the Spitzer source.The spectral index of the quiescent component is α=0.37±0.34between 6cm and 3.6cm.A factor of two increase in 6cm emission was detected during one epoch and circular polarization was marginally detected at the 5σlevel with Stokes V/I =48±16%.We have searched for 22GHz H 2O maser emission toward L1014-IRS,but no masers were detected during 7epochs of observations between June 2004and December 2006.L1014-IRS appears to be a low-mass,accreting protostar which exhibits cm emission from a thermal jet or a wind,with a variable non-thermal emission component.The quiescent cm radio emission is noticeably above the correlation of 3.6cm and 6cm luminosity versus bolometric luminosity,indicating more radio emission than expected.In this paper,we characterize the cm continuum emission in terms of observations of other low-mass protostars,including updated correlations of centimeter continuum emission with bolometric luminosity and outﬂow force,and discuss the implications of recent larger distance estimates on the physical attributes of the protostar and dense molecular core.Subject headings:radiation mechanisms:thermal,non-thermal —radio continuum:stars —stars:for-mation1.INTRODUCTIONIt is extremely diﬃcult to identify the incipient stages of low-mass (≈1M sun )star formation be-cause dense molecular cloud cores obscure nascent pro-tostars.Submillimeter dust continuum surveys (e.g.,Ward-Thompson et al.1994,Shirley et al.2000,Visser et al.2002,Kirk et al.2005)have identiﬁed several dense cores with no apparent internal sources,based on the lack of an IRAS point source and the diﬀuse na-ture of submillimeter dust emission.Observations by the Cores to Disk Legacy Team (c2d)with the Spitzer Space Telescope have identiﬁed a few mid-infrared sources that are embedded near the submillimeter continuum peaks of previously classiﬁed starless cores (e.g.Young et al.2004,Bourke et al.2006).These new objects are of low-luminosity (L int ≤0.1L ⊙)and presumably low-mass since they were not previously detected by IRAS.Some of these objects may be in the earliest stages of accre-tion.These newly identiﬁed low-mass,low-luminosity protostars warrant detailed follow-up studies to deter-mine their evolutionary status.The ﬁrst newly identiﬁed object detected in the c2d survey,L1014-IRS (Young,et al.2004),was mod-eled as a very low-luminosity (L int <0.1L ⊙),low-mass1Bart J.Bok Fellow,Steward Observatory,University of Ari-zona,933Cherry Ave.,Tucson,AZ 857212NRAO,P.O.Box 0,1003Lopezville Road,Socorro,NM 878013Harvard-Smithsonian Center for Astrophysics,60Garden St.MS42,Cambridge,MA 021384Nicholls State University,Thibodaux,LA 703105Division of Geological and Planetary Sciences 150-21,Califor-nia Institute of Technology,Pasadena,CA 91125(M <0.1M ⊙)object embedded within the Lynds 1014dark cloud (Lynds 1962)at a distance of approximately 200pc.This object has been classiﬁed as a VeLLO (Very Low-Luminosity Object)by the c2d team:an object with an internal protostellar luminosity ≤0.1L ⊙that is di-rectly associated with a dense molecular core.The recent study of Morita et al.(2006)has suggested a revised dis-tance estimate of 400to 900pc based on the possible age ranges of nearby T-Tauri stars that are spatially within 2◦of the L1014dense core.However,it is not clear that these T-Tauri stars are directly associated with L1014.Determining the evolutionary state of L1014-IRS has been the subject of several follow-up studies.The large scale molecular distribution in the dense core was deter-mined by the single-dish mapping survey of Crapsi et al.(2005).No evidence for a large scale CO outﬂow was detected;however,a molecular outﬂow was detected on small scales with the SMA (Bourke et al.2005).Near-infrared observations detect scattered light,presumably from the outﬂow cone,at 1.6µm and 2.2µm (Huard et al.2006).The SMA-detected CO outﬂow is aligned with the direction of the near-infrared scattered light cavity.High resolution molecular observations with BIMA indi-cate that the protostar is not at the peak of the molec-ular and dust column density in the core,but oﬀset by about 8′′in the plane of the sky (Lai et al.,in prepa-ration).This oﬀset is also seen in (sub)millimeter con-tinuum maps (Young et al.2004)and the near-infrared extinction map (Huard et al.2006).Despite these signiﬁcant observational eﬀorts,a single,consistent picture of the evolutionary state of L1014-IRS has not emerged.In order to better characterize the2physical nature of L1014-IRS,we have conducted cen-timeter radio continuum observations using ﬁve array conﬁgurations of the Very Large Array 6(VLA).Cen-timeter radio continuum emission is well correlated with the luminosity of protostellar sources (Anglada 1995)and is thought to arise from shock ionization from protostel-lar winds (Ghavamian &Hartigan 1998),from interac-tion of the protostellar jets with dense gas in the interface of the outﬂow cavity (Curiel et al.1987,1989,Shang et al.2004),or from accretion shock-driven photoionization (Neufeld &Hollenbach 1996).In this paper we report the detection and characteriza-tion of the cm radio continuum emission toward L1014-IRS (§3.1).We compare the detected cm emission with observations of other low-mass protostars (§4.1).We dis-cuss the non-detections of 22GHz H 2O masers and com-pare our upper limits to the recent maser monitoring surveys of low-mass protostars (§4.2).We also compare the derived source properties from the diverse studies of the protostar and dense core in terms of the range of distance estimates to L1014(§4.3).2.VLA OBSERVATIONSL1014-IRS was observed during 10epochs in 5ar-ray conﬁgurations (D,A,BnA,B,and C)with the Very Large Array (Table 1).All observations were centered on the published Spitzer mid-infrared source (α=21h 24m 07s .51,δ=+49◦59′09.′′0,J2000.0).Con-tinuum observations were made at 3.6cm,and 6.0cm,with two polarization pairs at adjacent frequencies,pro-viding a total equivalent bandwidth of 172MHz.We also attempted to detect H 2O masers by observing the J K a K c =616→523transition at 22.23508GHz with,typically,24.4kHz spectral resolution (0.3km/s)span-ning ±20km/s velocity coverage.For the 3.6and 6.0cm data,the data were reduced independently using the standard routines in AIPS++and AIPS .Complex gain calibration was performed by switching to the nearby quasar 2137+510,2.4◦from L1014-IRS,on time-scales of 15to 30minutes (Hamaker,Bregman,&Sault 1996a,b).The absolute ﬂux density and bandpass calibration were determined from observa-tions of the quasars 3C48and 3C283.The Stokes I and V images were deconvolved using the Cotton-Schwab algorithm (e.g.,Schwab 1984)and Clark-H¨o gbom algorithm (H¨o gbom 1974,Clark 1980)with a few thousand iterations and interactive CLEAN regions.Imaging the L1014-IRS ﬁeld was diﬃcult due to the presence of several bright sources within the VLA pri-mary beam (Figure 1a).Special care had to be taken in the CLEANing process (e.g.,Cornwell,Braun,&Briggs 1999)and multiple reductions with variations in the CLEAN parameters were performed.We have checked the consistency of our images by also reducing the data with the standard AIPS routines,and the ﬂuxes agree within the statistical errorbars.Generally,the images are made with natural weighting of the visibilities (Briggs 1995);however,uniform weighting was used to obtain better angular resolution for the full track (9hour)ob-servations on the days of July 1,2004and August 21,6The VLA is operated by NRAO.The National Radio Astron-omy Observatory is a facility of the National Science Foundation operated under a cooperative agreement by Associated Universi-ties,Inc.2004.3.RESULTS3.1.Radio Continuum DetectionsWe detected cm radio continuum emission from a source within 0.′′2of the Spitzer mid-IR source using the VLA at 3.6cm and 6cm (Figure 1).The initial detec-tions were made during 9hour tracks at 3.6cm and 6cm with the VLA in the D conﬁguration.Subsequent observations detected the source at 6cm in the other three VLA conﬁgurations,and again at 3.6cm in D con-ﬁguration.All of the 6cm observations,except for the initial 6cm detection on August 21,2004,indicate a constant ﬂux density source with an average 6cm ﬂux of 88±11µJy (8σ).The source is unresolved in all array conﬁgurations.The two detections at 3.6cm are also in agreement with each other despite being separated by 17months.The average 3.6cm ﬂux is 111±8µJy (14σ).The spectral index between two wavelengths (λ2>λ1)is deﬁned asα=ln(S λ1/S λ2)3to detect that variation due to the higher noise level in progressively shorter time blocks.A Stokes V source (circularly polarized)was marginally detected at the 5σlevel at 6cm on August 21,2004with a ﬂux of 84±17µJy at the position of the 6cm Stokes I source (Figure 1d).The fraction of circular polarization,f c =Stokes V/I =48±16%,is quite high if the detection is signiﬁcant.Since the Stokes V result was obtained during a full single track of VLA observa-tions,it will be diﬃcult to conﬁrm until L1014-IRS is observed with wider bandwidth (e.g.,with the eVLA).If the circular polarization signal is real,then this observa-tion indicates that the radio emission observed on August 21,2004must originate from a non-thermal mechanism (§4.1.2).3.2.Water Maser SearchWe searched for H 2O masers during seven epochs span-ning 22months by observing the J K a K c =616→523transition at 22.23508GHz.No H 2O masers were de-tected at any epoch.The combined 1σrms of the non-detection is 3.1mJy/beam with a channel spacing of 24.4kHz and a total bandwidth of 3.125MHz (∼40km/s).The individual observations are summarized in Table 1.4.DISCUSSION4.1.Centimeter Radio Continuum Emission4.1.1.Steady ComponentCentimeter continuum emission has been detected to-ward many but not all high-mass and low-mass proto-stars.The emission is most commonly thought to origi-nate from bremsstrahlung (free-free)emission from ion-ized gas,although some protostellar objects also display non-thermal emission.For high-mass protostars,the ion-ization mechanism is photoionization usually in the form of an embedded HII region (Churchwell 1990).For pro-tostars that are later in spectral type than B,the ioniz-ing radiation from the star is not enough to signiﬁcantly photoionize the surrounding envelope and an alternative mechanism is needed to explain the observed emission (Rodr´ıguez et al.1989,Anglada 1995).Since nearly all low-mass,embedded protostars are known to have molecular outﬂows (e.g.Wu et al.2004),the ioniza-tion is postulated to arise from shocks generated from a jet (e.g.,Cohen,Bieging,&Schwartz 1982,Bieging &Cohen 1989,Curiel et al.1987,1989,Rodr´ıguez &Reipurth 1996,Shang et al.2004).Direct evidence for this hypothesis comes from observations using interfer-ometers at high angular resolution.Elongated centime-ter continuum emission regions are observed with the same orientation as the large-scale molecular outﬂow to-ward a few protostars (e.g.,Anglada 1995;Bontemps,Ward-Thompson,&Andr´e 1996).In addition,the ob-served outﬂow force (M ⊙km/s/yr)theoretically provides enough energy in the shock to explain the observed cen-timeter ﬂuxes toward most low-mass protostars (Cabrit &Bertout 1992,Skinner et al.1993,Anglada 1995).Finally,the radio spectral index of many protostellar sources is consistent with optically thin (α≈−0.1)to partially optically thick free-free emission between 3.6cm and 6.0cm (2.0>α>−0.1,e.g.,Anglada et al.1998,Beltr´a n et al.2001).The centimeter continuum luminosity (e.g.,L 3.6=S 3.6D 2mJy kpc 2)of low-mass and intermediate-mass protostellar sources was ﬁrst cataloged from the litera-ture in the review by Anglada (1995).Anglada plotted the 3.6cm luminosity against the bolometric luminosity of protostars with L bol <103L ⊙and found a well corre-lated relationship (r =0.79),L 3.6=10−2.1(L bol /1L ⊙)0.7mJy kpc 2.This relationship has formed the basis for predictions of the amount of centimeter emission that is expected in searches for new low-mass protostars (e.g.,Harvey et al.2002,Stamatellos et al.2007).The 3.6cm luminosity correlation is directly related to the well established correlation of outﬂow force vs.bolometric lu-minosity (Bontemps et al.1996,Wu et al.2004);higher luminosity sources drive more powerful outﬂows that re-sult in a larger degree of shock ionization and therefore a larger 3.6cm luminosity (Anglada 1995).Since 1995,many more centimeter observations have been made and the spectral coverage of the photometry of protostars has increased.We have used the detailed summary tables of Furuya et al.(2003;Table 4)and Anglada (1995)supplemented by the surveys of Eiroa et al.(2005)and Anglada et al.(1998)to catalog the 3.6cm,6.0cm,and bolometric luminosities of detected protostellar sources.We updated the bolometric lumi-nosity of sources observed in the submillimeter surveys of Shirley et al.(2000),Mueller et al.(2003),and Young et al (2003).The resulting sample of 58sources at 3.6cm and 40sources at 6.0cm are plotted in Figure 2.We ﬁnd updated correlations oflog(L 3.6/1mJy kpc 2)=−(2.24±0.03)+(0.71±0.01)log(L bol /1L ⊙)(2)log(L 6.0/1mJy kpc 2)=−(2.51±0.03)+(0.87±0.02)log(L bol /1L ⊙),(3)with correlation coeﬃcients of r =0.66and r =0.74respectively.This sample is not complete as there are many more protostellar centimeter detections for which no L bol has been published.Nevertheless,we have up-dated the correlation of Anglada (1995)with twice as many points at 3.6cm and plotted the correlation at 6.0cm for the ﬁrst time.For comparison,we have plotted the 3.6cm and 6.0cm luminosities of L1014-IRS in Figure 2at the distances of 200,400,and 900pc.L1014-IRS is above the correlation at all distances indicating that we have detected more centimeter continuum ﬂux than expected.The spectral index of sources between 3.6cm and 6.0cm is used to argue for the interpretation that the emission mechanism is consistent with partially optically thick free-free emission.Optically thin free-free emission is expected to have a spectral index α=−0.1at cen-timeter wavelengths.In the optically thick limit,αap-proaches 2.0with intermediate values indicative of par-tially optically thick plasmas.We have also plotted the spectral index of sources from the literature that have been detected at both wavelengths and have a published L bol determination in Figure 2.No correlation of αis observed with bolometric luminosity,probably indicat-ing that the optical depth associated with the jet’s shock ionization is not dependent on the total protostellar lu-minosity or the strength of the molecular outﬂow.The median spectral index is α=0.5,with most protostellar sources having ﬂat or positive spectral indices.This me-dian value is close to the result expected for an ionized wind or jet with a 1/r 2density gradient (e.g.,Panagia &Felli 1975,Wright &Barlow 1975,Reynolds 1986).Unfortunately,most of the sources in Figure 2were not4observed at both wavelengths on the same day and vari-ability may result in signiﬁcant scatter in the plot.The spectral index of the quiescent emission of L1014-IRS(0.37±0.34)is consistent with ionized free-free emis-sion with a density gradient.The spectral index agrees well with the median of the ensemble of protostellar sources measured.We have also updated the correlation between outﬂow force,F out(M⊙km/s/yr),and centrimetric luminosity of Anglada(1995)using the molecular outﬂow compilations of Bontempts et al.(1996),Furuya et al.(2003),and Wu et al.(2004).The updated correlation of44sources is weak(r=0.55),log(F out/1M⊙km/s/yr)=−(3.15±0.07)+(0.67±0.03)log(L3.6/1mJy kpc2).(4) This correlation is interpreted as evidence that jets from the molecular outﬂow provide enough shock ionization to account to the observed centimeter continuum emission (e.g.,Anglada1995).Assuming maximum ionization ef-ﬁciency,the minimum outﬂow force needed to create the observed level of3.6cmﬂux was estimated by Curiel et al.(1987,1989)to be F out=10−3.5(L3.6/1mJy kpc2). This level is shown as a dashed line in Figure2d.The observed outﬂow force toward L1014-IRS ranges from 0.04−2.9×10−6M⊙km/s/yr for distances of200to 900pc(see Bourke et al.2005).At a distance of200pc, the upper limit on the outﬂow force is a factor of2lower than the Curiel theoretical minimum outﬂow force.The upper limit in the observed outﬂow force includes esti-mates of the missingﬂux due to interferometric spatial ﬁltering as well as the average opacity corrections for low-mass protostellar outﬂows(see Bourke et al.2005). Given the range of uncertainty in these estimates,the observed outﬂow force is below,but not necessarily in-consistent with the theoretical minimum outﬂow force needed to produce the observed3.6cm luminosity.The disagreement between the observed outﬂow force is more pronounced at larger distances,increasing to an order of magnitude below the theoretical curve for a distance of 900pc.The small observed molecular outﬂow force and the large observed centimeter continuum luminosities may indicate that another ionization mechanism is operating in L1014-IRS.There are a few other possibilities that have been discussed in the literature.We shall analyze the viability of two popular possibilities.Radiative transfer modeling of emission at24and70µm indicate aﬂux excess from a disk around L1014-IRS (Young et al.2004).Neufeld&Hollenbach(1996)pos-tulated that the supersonic infall of material onto a pro-tostellar disk will create an accretion shock with enough ionization to generate∼1mJy of continuum emission at centimeter wavelengths at distances of≈200pc.How-ever,this mechanism does not appear to be able to pro-vide enough ionization to explain the emission observed toward L1014-IRS.In order to produce90µJy emission at6cm at a distance of200pc,the protostellar mass has to be>2M⊙and the accretion rate onto the disk must be>10−4M⊙/yr(see Figure2of Neufeld&Hollenbach 1996).Estimates of the protostellar mass are very un-certain and highly distance dependent,however,2M⊙is likely larger than the L1014-IRS protostar and disk mass, even for the far distance estimate of900pc(see Young et al.2004).Furthermore,this accretion rate is an order of magnitude larger than the range of inferred accretion rates from the observed outﬂow momentumﬂux and the modeled internal luminosity of the source(≤3×10−5 M⊙/yr;Bourke et al.2006,Young et al.2004).Ion-ization from an accretion shock does not appear to be a likely explanation.A more plausible possibility is that there is a spheri-cal wind component.The expected continuum emission from self-shocked spherical winds have been modeled by numerous authors(e.g.,Panagia&Felli1975,Wright& Barlow1975,Reynolds1986,Gonz´a lez&Cant´o2002). The recent study of Gonz´a lez&Cant´o model a time variable wind that generates internal shocks(Raga et al. 1990)which produce ionization and centimeter contin-uum emission(see Ghavamian&Hartigan1998).Their models produce centimeter continuum emission of sev-eral hundreds ofµ-Janskys and spectral indices that are positive for mass loss rates of10−6M⊙/yr at a distance of150pc.Accounting for the larger distance estimates of200to900pc and potentially lower mass loss rates for L1014-IRS,then this type of emission may still account for theﬂux observed at3.6cm and6cm(∼100µJy). However,the spherical wind models is usually applied to evolved protostellar sources that are Class II(classi-cal T-Tauri stars)or later(e.g.,Evans et al.1987).If L1014-IRS is an older,more evolved protostar,then this may not be a problem(§4.3).4.1.2.Variable ComponentWhile the quiescent component of L1014-IRS is constant over5epochs with a positive spectral in-dex,the emission properties of L1014-IRS were signif-icantly diﬀerent during the single epoch of August21, 2004.The6cmﬂux was larger by a factor of two, S6.0(21AUG2004)=173±16µJy.Unfortunately,the spectral index of elevated emission was not determined since L1014-IRS was observed at a single wavelength. Circular polarization was detected at the5σlevel indi-cating non-thermal emission.Thus,L1014-IRS has vari-able centimeter emission,although the timescale of the variability is not constrained since it was seen to vary during only a single epoch.In general variability of centimeter continuum sources has not been properly addressed since observations of sources are limited to a few epochs.There has not been a systematic,monthly monitoring campaign of deeply embedded sources to characterize their centimeter vari-ability;however,there is observational evidence for vari-ability among deeply embedded protostellar sources.For instance,the Class0source,B335,is known to vary be-tween an upper limit of<80µJy(1994December)and 390µJy(2001January)at3.6cm(Avila et al.2001, Reipurth et al.2002).Another example is the variability and purported evidence for jet precession of the centime-ter continuum sources toward IRAS16293(Chandler et al.2005).Variability contributes scatter(de-correlation) of the luminosity correlations shown in Figure2.In ad-ditional to the known variable thermal sources,several non-thermal protostellar sources are known to be highly variable(e.g.T Tauri stellarﬂares,see White1996); but,most of those objects are more evolved than deeply embedded protostars.The observed increase in emission on August21,2004 and5σStokes V detection is indicative of variable non-5thermal emission toward L1014-IRS.While there have been a few high-mass protostars with observed nega-tive spectral indices(e.g.,Reid et al.1995;Garay et al.1996),most embedded(≤Class I)low-mass proto-stars with cm radio emission have positive spectral in-dices(see Figure2).There are only a few embedded low-mass protostars toward which negative spectral in-dex,non-thermal emission is detected(e.g.,Shepherd &Kurtz1999,Girart et al.2002).One case,R CrA IRS5,was detected with signiﬁcant circular polarization (Feigelson et al.1998).The authors postulate that the emission is due to gyrosynchrotron emission and may originate from magnetic reconnection events associated withﬂares(Feigelson et al.1998).The physical mecha-nism for generating the radioﬂare is still not well under-stood(e.g.,Basri2004)and it is questionable whether it applies to the embedded phase of low-mass protostars (see below).A second case,IRAS19243+2350,is a steep spectrum non-thermal source(α=−0.82±0.04)that is elongated in the direction of its CO outﬂow(Girart et al.2002).Girart et al.postulate that the emission may originate from a bi-conical synchrotron source tracing the protostellar jet,similar to observations of the high-mass source W3(OH)(Reid et al.1995,Wilner,Reid,& Menten1999).Non-thermal emission may be present in low-mass protostellar jets,but the level of emission may be dominated by the thermal,shock-ionized component of the jet.In the case of L1014-IRS,this non-thermal jet component would have to be variable.We shall dis-cuss several possibilities for the origin of the observed non-thermal emission toward L1014-IRS.In order to better understand the non-thermal emis-sion mechanism,we estimate the brightness tempera-ture of the emission to be T b=Sνλ2d2/2kR2emit= 0.14−2.8×106K for distances of200to900pc and an emitting region that is R emit=1AU in size.The brightness temperature is very sensitive to the assumed size of the emitting region.A R emit of1AU is appro-priate for a smallﬂare;but,smaller when compared to the solar coronal emitting region which is typically less than≈5AU(e.g.Leto et al.2000).Since the emission observed toward L1014-IRS was unresolved,even in the A-array conﬁguration,then we can only limit the size of the emission region to<90(D/200pc)AU.For instance, if we assumed that the emitting region was45AU(half of our A-array resolution),then the brightness temper-ature drops to<100K.This is too low even for the steady,thermal free-free component(T∼104K),unless the emission was very optically thick.This cannot be the case sinceτ>>1would implyαapproaching2.0 which is not observed.While the size of the emitting re-gion is severely unconstrained,the detection of circular polarization indicates non-thermal emission probably on small(<few AU)size scales,most likely due to gyrosyn-chrotron emission(e.g.,Ramaty1969;Dulk&Marsh 1982).Radioﬂares are observed toward very low-mass ob-jects including brown dwarfs(Berger et al.2002,Os-ten et al.2006),late M dwarfs(Berger2006),and T-Tauri stars(Bieging&Cohen1989;White,Pallavicini, &Kundu1992).The typical level of quiescent emis-sion toward low-mass stars and brown dwarfs is100µJy and theﬂares are1mJy(for nearby distances of≈30pc)with signiﬁcant circular polarization(f c≈50%) detected during theﬂaring events.The spectral lumi-nosity at6cm of L1014-IRS during the event is Lν= 4πD2Sν=4.7×1015(D/200pc)2erg s−1Hz−1.This is about2400times larger than the most luminous ob-served solarﬂares(Bastian2004)and about200times larger than theﬂares observed by Berger(2006)toward late M dwarfs.If the radio emission is due to aﬂare,it must be a powerfulﬂare since L1014-IRS is at least20 times farther away than the average distance of sources detected by Berger( d =10.6±5.1pc).However,it is not larger than the typical non-thermalﬂaring emission observed toward young T Tauri stars of Lν=1015−1018 erg s−1Hz−1(G¨u del2002).The ratio of spectralﬂare luminosity to bolometric luminosity,Lν/L bol=4×10−18 Hz−1,is also similar to the ratio observed toward classi-cal T-Tauri stars(see G¨u del2002).The timescale over which radioﬂares toward low-mass stars are observed tends to occur over minutes to hours. For instance,the low-mass brown dwarf,LP944-2,dis-covered by the2001NRAO summer students(Berger et al.2002),displaysﬂaring activity with an average timescale of10to15minutes.This is very diﬀerent from the activity observed toward L1014-IRS on August21, 2004.We detected no evidence for short-term variabil-ity within our detected emission.The source appears to have a nearly constantﬂux that is twice as high as the steady component for at least an8hour period.This elevated emission then appears to be longer in duration than that observed duringﬂaring events toward low-mass (proto)stars(G¨u del2002).An alternative possibility is that the elevated emission is not due to aﬂare,but due to rotational modulation of a non-thermal component associated with the magnetic connection between disk and accretion onto the star(i.e., Bieging&Cohen1989).Such a mechanism has been postulated for the T-Tauri star,V410Tauri,with a rota-tional period of1.9days.The observed emission toward V410is1mJy with a negative spectral index.If the accretion spot is blocked from view for a fraction of the stellar rotation period,then it is possibly that we could have observed the source with the accretion spot is in view on August21,2004and with the accretion spot blocked from view during the other epochs.The rota-tional period of L1014-IRS must be longer than8hours since elevated emission was observed during the entire8 hour track.A negative spectral index was observed to-ward V410Tauri,while a negative spectral has not been observed toward L1014-IRS.This hypothesis is highly speculative and would require a regular monitoring cam-paign to test.Unfortunately,it is currently not possible to strongly constrain the origin of the non-thermal component to-ward L1014-IRS.The expanded bandwidth of the eVLA is needed to permit a systematic monitoring campaign with high enough signal-to-noise in only a few hour ob-servations to routinely check for a Stokes V detection and to determine an instantaneous spectral index.4.2.Water Maser Non-detectionsA compact,weak molecular outﬂow has been detected toward L1014-IRS(Bourke et al.2005);therefore it may be possible to detect water masers if the jet is impinging on dense knots of material near the protostar.Previous。

Robust general equilibrium under stochastic volatility model

dDt ¼ lDt dt þ dv t
pﬃﬃﬃﬃﬃ v t Dt dBt ; pﬃﬃﬃﬃﬃpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ jðh À v t Þdt þ rv v t 1 À q2 dZ t þ qdBt ;
ð1Þ ð2Þ
where Bt and Zt are standard independent Brownian motions. This endowment ﬂow model is the standard stochastic volatility model of Heston (1993) with constant mean growth rate l P 0 and local variance vt.1 vt is a square-root mean reverting process with long-run mean h, speed of adjustment j, and variation coefﬁcient rv of the diffusion volatility vt. The parameters j,h, and rv are assumed to be nonnegative and satisfy 2jh P r2 v. Fixing the time period at D, the agent is assumed to maximize the following intertemporally additive expected utility (in the absence of a preference for robustness)
W. Xu et al. / Finance Research Letters 7 (2010) 224–231

Domsch2013_Increased BOLD sensitivity in the orbitofrontal cortex

Increased BOLD sensitivity in the orbitofrontal cortex usingslice-dependent echo times at 3T ☆Sebastian Domsch a,⁎,1,Julia Linke b,1,Patrick M.Heiler a ,Alexander Kroll a,b ,Herta Flor b ,Michèle Wessa b ,Lothar R.Schad aaDepartment of Computer Assisted Clinical Medicine,Medical Faculty Mannheim,Heidelberg University,Theodor-Kutzer-Ufer 1–3,68167Mannheim,GermanybDepartment of Cognitive and Clinical Neuroscience,Central Institute of Mental Health,Medical Faculty Mannheim,Heidelberg University,GermanyReceived 13December 2011;revised 12June 2012;accepted 21June 2012AbstractFunctional magnetic resonance imaging (fMRI)exploits the blood oxygenation level dependent (BOLD)effect to detect neuronal activation related to various experimental paradigms.Some of these,such as reversal learning,involve the orbitofrontal cortex and its interaction with other brain regions like the amygdala,striatum or dorsolateral prefrontal cortex.These paradigms are commonly investigated with event-related methods and gradient echo-planar imaging (EPI)with short echo time of 27ms.However,susceptibility-induced signal losses and image distortions in the orbitofrontal cortex are still a problem for this optimized sequence as this brain region consists of several slices with different optimal echo times.An EPI sequence with slice-dependent echo times is suitable to maximize BOLD sensitivity in all slices and might thus improve signal detection in the orbitofrontal cortex.To test this hypothesis,we first optimized echo times via BOLD sensitivity simulation.Second,we measured 12healthy volunteers using a standard EPI sequence with an echo time of 27ms and a modified EPI sequence with echo times ranging from 22ms to 47ms.In the orbitofrontal cortex,the number of activated voxels increased from 87±44to 549±83and the maximal t -value increased from 4.4±0.3to 5.4±0.3when the modified EPI was used.We conclude that an EPI with slice-dependent echo times may be a valuable tool to mitigate susceptibility artifacts in event-related whole-brain fMRI studies with a focus on the orbitofrontal cortex.©2013Elsevier Inc.All rights reserved.Keywords:Echo time;Event-related fMRI;Susceptibility gradients;BOLD sensitivity;Orbitofrontal cortex1.IntroductionIn the past two decades,functional magnetic resonance imaging (fMRI)exploiting the blood oxygenation level dependent (BOLD)effect [1,2]has revealed valuable in-sights of the brain functions underlying psychological processes like decision making and reward processing [3].The BOLD effect associated with oxy-and deoxyhemoglo-bin and perfusion changes [4]has been analytically de-scribed [5,6]and tested in phantoms [7].Neuronal activationin orbitofrontal regions is dif ﬁcult to image as signal loss and image distortions are common in this brain region [8].These detrimental effects result from macroscopic suscepti-bility gradients occurring close to air/tissue interfaces and are a general problem of the gradient-echo echo-planar imaging (EPI)sequence [9]commonly used for fMRI [10,11].Several strategies to compensate for susceptibility effects have been proposed and include z -shimming [12–14],tailored excitation pulses [15–17]or spiral trajectories [18–20].Either these strategies require sophisticated mod-i ﬁcations of the pulse sequences and processing routines or they compromise spatial or temporal resolution.An alternative strategy is the optimization of echo times since the BOLD contrast depends critically on this sequence parameter particularly in subcortical brain areas [21].Robust activation in the amygdala,also affected by susceptibility gradients,and the orbitofrontal cortex has been observed inAvailable online at Magnetic Resonance Imaging 31(2013)201–211☆Con ﬂict of interest declaration:None of these authors has any current or potential con ﬂict of interest concerning this paper.⁎Corresponding author.Tel.:+496213835120;fax:+496213835123.E-mail address:sebastian.domsch@medma.uni-heidelberg.de (S.Domsch).1These authors contributed equally to this work.0730-725X/$–see front matter ©2013Elsevier Inc.All rights reserved./10.1016/j.mri.2012.06.020many previous whole-brain fMRI studies using an EPI sequence with a short echo time of 27ms [22–24].However,the orbitofrontal cortex extends over several slices,and each slice is affected differently by susceptibility gradients.For a block design focusing at the amygdala,Stoecker et al.[25]showed that an EPI sequence with slice-dependent echo times improved signal detection at 1.5T.This method leaves temporal resolution unaffected,which is important for event-related designs.Hence,it would be of great relevance for neuroimaging research if this approach could also improve signal detection in the orbitofrontal cortex during event-related fMRI compared to a standard EPI protocol with reduced echo time at ﬁeld strength of 3T.The goal of the present work was the improvement of BOLD sensitivity (BS)in the orbitofrontal cortex usingslice-dependent echo times.First,T 2*relaxation times and ﬁeld gradient measurements were performed to simulate BS and to approximate optimal echo times in single slices.Second,a standard EPI sequence with short echo time was compared to a modi ﬁed EPI sequence with slice-dependent echo times during an event-related whole-brain fMRI experiment with a focus on the orbitofrontal cortex.The practical relevance of the modi ﬁed EPI was judged by comparing the number of activated voxels and the statistical signi ﬁcance of the activation.2.Materials and methods 2.1.BS calculationDuring neuronal activation,the BOLD effect alters microscopic background ﬁelds and thus the relaxation rate T 2*.In order to detect local changes of neuronal activation,the signal has to be susceptible to small changes of T 2*.The BS is a measure of the local signal intensity change (dI)related to changes of local T 2*relaxation times dT 2*and is proportional to the derivative of the signal intensity withrespect to T 2*,BS ∝dI dT *2.Neglecting macroscopic back-ground ﬁelds,the signal intensity is given by I =ρ:exp −TE =T 2*ðÞ½ :ð1ÞThe ﬁrst-order Taylor expansion of the signal intensityevaluated at the T 2*value during the baseline state (T *2B )can be written as:I =ρ⋅e −TE =T *22B 1+TE T *22B ⋅T *2−T *2B !;ð2Þwhere ρis the local spin density.Eq.(2)motivates thefollowing de ﬁnition for the BS:BS :=T *22B dI dT *2T *2=T *2B:ð3Þ[21].Calculating the BS according to this de ﬁnition yields the simple expression:BS =TE *I ;ð4Þwhich has been experimentally veri ﬁed [26].Accordingly,the maximal BS is reached when TE is equal to T 2*,which is in accordance with literature results presented in the work of Gati et al.[27].In the following,the effects of macroscopic ﬁeld inhomogeneities on BS will be considered.Through-plane gradients (G ss )in the z -direction cause extra spin dephasing along the slice direction.Accordingly,the expression for the signal intensity given by Eq.(1)has to be modi ﬁed to:I =ρ:exp −TE =T 2*ðÞ½ ∫P z 0ðÞe i γG ss z 0TE dz 0;ð5Þwhere P(z)denotes the slice pro ﬁle and γthe gyromag-netic ratio.Assuming rectangular slice pro ﬁles,which were used for all measurements in this study,the integral in Eq.(5)can be analytically computed and is given by sinc [γ·G ss ·TE·Δz /2],where Δz denotes the slice thickness.As discussed in the work of Deichmann et al.[21],macroscopic in-plane gradients (G sp )parallel to the EPI phase encode direction alter the sampling rate and therefore also the signal intensity.They further shift the effective echo time (TE eff )of the central gradient echo of the EPI readout.Both effects can be described by the single variable Q ,de ﬁned as:Q =1−γΔt FoV G sp =2πðÞ;ð6Þwhere Δt denotes the time between consecutive gradient echoes and FOV denotes the ﬁeld of view in the EPI phase encoding direction.Q affects the signal and the echo time in the following way:I =I 0=Q and TE eff =TE =Q :ð7ÞIn Eq.(7),I 0is the signal intensity and TE is the effectiveecho time without in-plane gradients.Eq.(6)shows that Q is equal to one if the ﬁeld gradients G sp are zero.Depending on the sign of the gradients,Q becomes smaller or greater than one.According to Eq.(7),the gradient echo is shifted to smaller or to higher TE values if the in-plane gradients are parallel or antiparallel to the EPI phase encoding direction.The signal and therefore also the BS are completely lost if the gradient echo is shifted outside the EPI readout window of the duration TA:j TE −TE =Q j N TA =2ð8Þ[25],where TA depends on the number of phase encoding steps and the bandwidth per pixel.According to Eq.(5)and Eq.(7),in the presence of macroscopic ﬁeld inhomogeneities202S.Domsch et al./Magnetic Resonance Imaging 31(2013)201–211along the EPI phase encoding direction and along the slice direction,the signal intensity is given by:I =ρ=Q :exp −TE =QT 2*ðÞ½ ∫P z 0ðÞe i γG ss z 0TE dz 0:ð9ÞApplying the BS de ﬁnition given in Eq.(3)to the ﬁrst-order Taylor expansion of Eq.(9)yields the following expression for the BS in the presence of macroscopic background ﬁelds:BS =TE =Q 2ρ:exp −TE =Q :T 2*ðÞ½:sinc γ:G ss :TE :Δz =2½ :ð10ÞSince BS is proportional to 1/Q 2,the maximal value of the BS increases for positive gradients (Q b 1)and decreases for negative gradients (Q N 1)(i.e.,in-plane gradients are parallel or antiparallel to the EPI phase encoding direction).In the case of nonzero in-plane gradients (G sp ≠0;G ss =0),the optimal echotime is shifted to TE=QT 2*.Nonzero through-plane gradients (G ss ≠0;G sp =0)shift the BS maximum to lower TE due to additional signal dephasing in slice direction.2.2.MR measurementsAll measurements were performed on a 3T whole-body MR system (Magnetom Trio,Siemens Healthcare,Erlangen,Germany)equipped with a 12-channel head coil.Written informed consent was obtained from all subjects who participated in this study.Gradient ﬁeld maps were acquired from six healthy subjects (27±1years,ﬁve males)using a multi-gradient echo 3D-FLASH sequence.After applying the RF-pulse,multiple gradient echoes were generated under the readout gradient slope spanning a range of echo times.The sequence was repeated until k-space was fully sampled.The slices were aligned along the anterior and posterior commis-sure (AC –PC)line and the following sequence parameters were used:40slices;matrix size=64×64;repetition time (TR)=30ms;ﬂip angle α=12°;FOV=220mm;slice thickness Δz =2.5mm.From each of the contrasts,one phase image was reconstructed,yielding a set of phase images at the different TE values of 2.45ms,5.12ms,7.79ms,10.46ms,and 13.13ms.The phase difference Δφbetween adjacent voxels depends on the strength of the macroscopic ﬁeld gradient G (assumed to be linear)and is proportional to the voxel dimension Δy and the time t between the signal excitation and the readout,Δφ=γ·G·Δy·t .Thus,the gradient strength could be determined from the set of phase images acquired at different echo times.Phase unwrapping was performed along the temporal dimension [28]before G ss and G sp were estimated by two separate linear least squares ﬁt on a voxel-by-voxel basis:Δφss =γ:G ss :Δz :t ;and ð10a ÞΔφsp =γ:G sp :Δy :t :ð10b ÞAfterwards,T 2*-weighted intensity images were acquired at 12different echo times between 5ms and 140ms usingthe multigradient echo 3D-FLASH sequence with the follow-ing sequence parameters:TE =[5,10,15,20,25,30,40,50,70,90,110,140]ms;α=25°.The remaining sequence parameters were kept constant.It has been shown thatexponential T 2*-ﬁtting without accounting for through-plane dephasing leads to an underestimation of T 2*[29].Therefore,G ss resulting from the phase mapping study was inserted into Eq.(5),which then was used as the ﬁtting function.Con ﬁdencebounds for local T 2*-ﬁtting and gradient ﬁeld mapping were setto 95%.Voxels with a relative ﬁt error of ΔT 2*/T 2*N 0.25andΔG /G N 0.25were not further processed.The resulting T 2*andgradient ﬁeld maps were coregistered to the T 2*-weighted intensity images,spatially normalized to standardized co-ordinates in the Montreal Neurological Institute space and ﬁnally resampled to a 96×96×40matrix size.2.3.fMRI experiment with focus on orbitofrontal cortex activationWithin an interval of 2weeks,12healthy participants (23±3years,six females)performed a self-paced probabi-listic reward-reversal task twice using either the modi ﬁed or the standard EPI sequence (40slices;descending slice direction;matrix size=96×96;TR=2.7s;bandwidth=1578Hz,ﬂip angle=90°;FOV=220mm;slice thickness=2.3mm,interslice gap=0.7mm,parallel imaging factor PI =2using GRAPPA [30].).Axial slices were aligned along the AC –PC.Both types of EPI sequences were counterbalanced to avoid ordering effects.The echo time of the standard EPI was set to 27ms.The echo time of the modi ﬁed EPI sequence increased linearly from TE=22ms to TE=47ms (Fig.1)in ascending slice direction (i.e.from slice 1to 27)by implementing a slice-dependent delay between the signal excitation and the EPI readout train.In the upper slices (i.e.from slice 28to 40),a ﬁxed delay of 25ms was imple-mented,yielding constant echo times of TE=47ms.Delaying the echo times increased the temporal resolution (i.e.TR for one EPI image stack)by the sum of the delays in all slices ∑slice ¼1 (40)Δτslice .Compared to the EPI with a ﬁxed TE of 27ms,the actual TR of the slice-dependent TE EPI is prolonged by approximately 0.6s (assuming the parameter settings from above).If TR is not below 2.5s,this temporal loss can be compensated by shortening the dead time between the EPI readout and the next RF excitation depending on the matrix size,TE and the bandwidth.Linearly increasing echo times were used to avoid abrupt image intensity changes in adjacent slices,which may otherwise corrupt image coregistration in the process of spatial normalization [25].The reward-reversal paradigm was identical to the paradigm used by Linke et al.[31].In brief,volunteers made a forced choice between two play cards on each trial.The “correct ”play card received an 8/2ratio of positive to negative feedback,whereas the “incorrect ”play card always received negative feedback.Feedback was provided in the form of winning (reward)or losing (punishment)a small amount of money.203S.Domsch et al./Magnetic Resonance Imaging 31(2013)201–211Imaging data were preprocessed and analyzed according to standard procedures using SPM8(http://www.ﬁ/spm ).First,functional scans were slice-timed and then realigned to the ﬁrst scan of each ing the modi ﬁed compared to the standard EPI sequence,the slices were not uniformly sampled in time.Therefore,the standard SPM slice timing routine was adjusted accounting for the actual acquisition timing of the modi ﬁed EPI sequence.During normalization,images were resampled every 2.3mm using sinc-interpolation and were smoothed with a 7×7×7-mm 3Gaussian kernel to decrease spatial noise.For the reward N pun-ishment condition,the BOLD signal changes were modeled ona voxel-by-voxel basis for every single subject (SPM ﬁrst-level analysis)according to the general linear model.3.Results 3.1.BS calculationFig.2shows T 2*-relaxation times below 10ms and susceptibility gradients of more than 200μTm −1in the lower part (slice 1)of the orbitofrontal cortex observed in a single subject.The upper part (slice 5)is clearly less affected bysusceptibility gradients leading to increasing T 2*and decreasing susceptibility gradients values.The results of all subjects are summarized in Table 1.In slice 1,T 2*was on average between 8±3ms and 31±7ms.In-(through-)plane gradients were in the range of −202±57μTm −1(−230±56μTm −1)and 288±63μTm −1(214±17μTm −1).In slice 5,T 2*was in the range of 19±3ms and 47±9ms.In-(through-)plane gradients were in the range of −135±61μTm −1(−209±36μTm −1)and 73±25μTm −1(63±25μTm −1).Fig.3shows the effect of typical gradient strengths on BSobserved in the orbitofrontal cortex.For example,given a T 2*value of 35ms and a susceptibility gradient of 150μTm −1in the slice (phase)encoding direction would shifts the BS peak from 35ms to 4ms (19ms).Assuming the BOLD signal is acquired at a TE value of 27ms,then only 87%(2%)of the maximal BS would remain.Fig.4depicts a simulation of the maximal TE values before total signal loss occurs due to in-plane gradients intheFig.1.Slice-dependent echo time data acquisition scheme.Slices are aligned along the anterior to the center of the posterior commissure (AC –PC).The orbitofrontal cortex is imaged with increasing echo times between 22ms and 37ms.Fig.2.T 2*-relaxation time and ﬁeld gradient maps in slice and phase encoding direction depicted for two different axial slices.The orbitofrontal cortex is marked by the ROI.The data were obtained from a single subject.204S.Domsch et al./Magnetic Resonance Imaging 31(2013)201–211range of [−350…350]μTm −1.According to Eq.(8),the maximal TE (TE max )was de ﬁned as TE max ≡1/2·TA/|1−1/Q|.It shows that a TE value of 22ms prevents total signal loss for any negative in-plane gradients in the range of [−350…0]μTm −1(Q N 1).For increasing positive gradient strengths in the range of [0…180]μTm −1(0b Q b 1),the maximal echo time decreases rapidly below 20ms.In the range of [180…350]μTm −1(Q b 0),TE max slightly increases but still remains below feasible TE values.Fig.5shows a BS calculation based on Eqs.(8)and (10)for typical in-and through-plane gradients.According to this simulation,a TE value of 20ms is adequate to recover BS forall possible combinations of through-plane gradients in the range of −350μTm −1to 350μTm −1and in-plane gradients in the range of −350μTm −1to the order of 100μTm −1.When TE is reduced from 40ms to 30ms and to 20ms,the percentage of dead pixels (i.e.BS /BS max b 0.1)decreases from 37%to 32%and to 8%in the left hemisphere of the G sp –G ss plane (i.e.,G sp b 0).This echo time reduction yields nearly no effect in the right hemisphere (i.e.,G sp N 0),decreasing the percentage of dead pixels from 44%to 43%and to 40%.Fig.6demonstrates that,in inferior slices,shortening TE effectively increases BS and reduces signal loss in the orbitofrontal cortex but decreases BS in areas not affected by susceptibility gradients.The average percentage of dead voxels in the orbitofrontal cortex decreases from 67%±15%to 39%±16%to 16%±4%in all subjects when TE is reduced from 40ms to 30ms to 20ms.A calculation of the average BS in the orbitofrontal cortex in all subjects reveals that BS is signi ﬁcantly enhanced when TE increases with the slice number (Fig.7).Optimal TE values simulated for all subjects are depicted in Fig.8.The BS is maximal for increasing TE values between 17±2ms (slice 1)and 39±2ms (slice 9)in foot-to-head direction.In areas not affected by susceptibility gradients,an optimal TE value of 46±3ms was measured,which is in good agreement with cortical values of about 40–50ms [32].3.2.FMRI experiment with focus on orbitofrontal cortex activationFig.9depicts two inferior slices acquired with the standard and the modi ﬁed EPI sequence.The orbitofrontal cortex clearly shows less signal loss when TE is reduced from 27ms to 22ms (Fig.9).The single-subject fMRI analysis shows orbitofrontal cortex activation in 5out of 12participants when the stan-Table 1Transverse relaxation times (T 2*)and susceptibility gradients in slice and phase encoding direction (G ss ,G sp )observed in two different slices of the orbitofrontal cortex Subject T 2*/ms G ss /μTm −1G sp /μTm −1a 1[8…32][−178…190][−124…209]2[10…25][−240…240][−250…186]3[6…29][−250…324][−255…219]4[5…32][−200…286][−214…212]5[11…44][−100…338][−274…239]6[5…23][−243…350][−261…216]Mean [8±3…31±7][−202±57…288±63][−230±56…214±17]b 1[17…58][−146…57][−34…90]2[20…41][−205…43][−154…79]3[15…35][−257…95][−180…65]4[22…55][−222…48][−104…27]5[20…47][−212…41][−204…81]6[18…47][−211…92][−135…94]Mean[19±3…47±9][−209±36…63±25][−135±61…73±25]In a and b,the minimal and maximal parameter values denote the 5and 95percentiles measured in slice 1and in slice 5,respectively.The standard deviation denotes the variation between the differentsubjects.Fig.3.Echo time dependence of the BS in the presence of typical susceptibility gradients occurring in the orbitofrontal cortex.The maximum is normalized to 1.205S.Domsch et al./Magnetic Resonance Imaging 31(2013)201–211dard EPI sequence was used (Table 2).With the modi ﬁed EPI sequence,activation was detected in all participants,the average number of activated voxels increased from 87±44to 549±83,and the maximal t value,averaged over all sub-jects,increased from 4.4±0.3to 5.4±0.3.At a signi ﬁcance level of α≤10−4(uncorrected),a second-level one-sample t -test revealed an average number of 55voxels to be activated in the orbitofrontal cortex using the standard EPI sequence (Fig.10).At the same threshold,553voxels were activated when the modi ﬁed EPI sequence was used.Activation in inferior and superior slices acquired at echo times below and,respectively,above 27ms was clearly more robust.Slices acquired at echo times of about 27ms showed similar activationin both EPI data sets.However,parietal activation is stronger at TE values above 27ms,showing superiority of the modi ﬁed EPI paring statistical inferences of both EPI sequences in a second-level paired t -test at a signi ﬁcance level of α≤10−3(uncorrected),stronger activation was detected in 303(0)voxels of the orbitofrontal cortex when the modi ﬁed (standard)EPI was used (Fig.11).4.DiscussionThe goal of the present work was to demonstrate that the adjustment of echo times in single EPI slices improves BSinFig.4.Simulation of maximum echo times before signal loss occurs due to macroscopic susceptibility gradients along the EPI phase encoding direction based on Eq.(8).As can be seen,the slope depends critically on the parameter Q de ﬁned in Eq.(7).Fig.5.BS maps in the presence of in-and through-plane gradients depicted for different echo times.The reference point at (0|0)is indicated by the white cross-hair.The BS maximum is normalized to 1and values below 10%are set to 0.This simulation is based on Eq.(8)and Eq.(10).The following sequenceparameters were assumed:T 2*=25ms,bandwidth=1578Hz,TA=30ms,FOV=2202mm 2,slice thickness=3mm.206S.Domsch et al./Magnetic Resonance Imaging 31(2013)201–211the orbitofrontal cortex.First,BS simulations were per-formed to approximate optimal echo times for single EPI slices.Second,an EPI sequence with slice-dependent echo times was compared to an optimized EPI sequence with constant short echo time during an event-related fMRI ex-periment with focus on the orbitofrontal cortex.Both simulations and functional results indicate that BS is increased when slice-adjusted echo times are used.In the inferior orbitofrontal cortex,gradient strengths more than 250μTm −1were observed,which are in agreement with the ﬁndings of DePan ﬁlis et al.[33].It was demonstrated that even moderate susceptibility gradi-ents in the order of 150μTm −1shift the maximal BS to echo times clearly below 30ms.Both BS simulations and actual EPI measurements con ﬁrmed that a TE value of about 20ms signi ﬁcantly reduces signal dropouts and increases BS in inferior slices of the orbitofrontal cortex.Further,weobserved that T 2*-relaxation times increase and ﬁeld gradients decrease in foot-to-head direction.Therefore,short TE values of about 20ms decrease the BS in superior slices of the orbitofrontal cortex.Accordingly,the numer-ically calculated optimal echo times averaged over all subjects increased with the slice number and were in the range of 17±2ms and 39±2ms.However,echo times below 20ms are not realizable using an EPI sequence with standard imaging parameters (e.g.matrix size,bandwidth,partial Fourier,parallel imaging factor,etc.).But using a short TE value of 20ms potentially preserves BS in the presence of strong in-and through-plane gradients of more than 300μTm −1.An important case-by-case analysis must be made with respect to the polarity of the in-plane gradients as described by Deichmann et al.[21,34].Our simulations showed that a TE value of 22ms is adequate to avoid signal loss caused by negative in-plane gradients of more than 300μTm −1,whereas an echo time reduction is ineffective to recover signal loss due to positive in-plane gradients.Even positive in-plane gradients of only 80μTm −1cannot be compensated via TE reduction because the gradient echo is rapidly shifted outside the acquisition window with increasing gradient strengths.We showed that an EPI with slice-dependent echo times between 22ms and 37ms signi ﬁcantly improves statistical inference in the orbitofrontal cortex in both on thesingle-Fig.6.BS maps obtained from a single subject and depicted for different echo times.The orbitofrontal cortex is marked by theROI.Fig.7.Average BS in ascending axial slices of the orbitofrontal cortex depicted for different echo times.The error bars denote the standard error of the mean.207S.Domsch et al./Magnetic Resonance Imaging 31(2013)201–211subject level and on the group level compared to an optimized standard EPI with short TE of 27ms.The effectiveness of the modi ﬁed EPI sequence might be further improved by changing the slice tilt to reduce in-plane gradients with positive signs as we have observed in the orbitofrontal cortex.It has been previously shown that a tiltangle of 30°from the axial towards the coronal orientation yields the best choice to minimize signal dropouts and thus to increase the BS in the orbitofrontal cortex [34].How-ever,it should be noted that,due to the shape of the brain,a transversal slice orientation yields time-ef ﬁcient volume coverage.The number of slices would potentially have to be increased when tilting the slice orientation away from a transversal toward a coronal orientation [34],which would be detrimental for the temporal resolution,which is particu-larly important in event-related fMRI studies.It has been experimentally demonstrated that the phase encoding scheme signi ﬁcantly affects susceptibility-induced signal losses in EPI [33].Hence,a slice-dependent change oftheFig.8.Optimal echo times in ascending axial slices of the orbitofrontal cortex and depicted for six healthyparticipants.Fig.9.Representative axial slices acquired at TE of 27ms (top)and about 22ms (bottom).Signal dropout in the orbitofrontal cortex is marked by the arrow.Table 2Single-subject fMRI results shown for the standard and the modi ﬁed EPI sequence SubjectStandard EPI Modi ﬁed EPI NT max N T max 1201 3.8289 4.42108 4.4359 4.83526 6.3609 6.040-461 5.1567 3.6127 3.960-1054 6.570-377 4.680-695 4.9977 5.0630 5.9100-243 4.81150 3.4873 6.4120-8667.3Mean87±444.4±0.3549±835.4±0.3N and T max denote the number of signi ﬁcant voxels and maximal t -values detected in the orbitofrontal cortex.The results are displayed at a threshold of α≤5·10−3(uncorrected)and a minimal cluster size of 15voxels.The error denotes the standard error of the mean.208S.Domsch et al./Magnetic Resonance Imaging 31(2013)201–211。

Volatility, Correlation and Tails for Systemic Risk Measurement

∗
1
Electronic copy available at: /1611229
1
Introduction
The Great Recession of 2007/2009 has motivated market participants, academics and regulators to better understand systemic risk. A useful deﬁnition of systemic risk by Federal Reserve Governor Daniel Tarullo, is “Financial institutions are systemically important if the failure of the ﬁrm to meet its obligations to creditors and customers would have signiﬁcant adverse consequences for the ﬁnancial system and the broader economy.” In this deﬁnition, it is the failure of a systemically important ﬁrm to meet obligations that is the cause of systemic distress. Thus, measures of systemic risk are associated with ﬁrm bankruptcies or near bankruptcies which are inevitable consequences of a decline in equity valuations for highly levered ﬁrms. This idea is further developed in the theoretical analysis of Acharya et al. (2010). A ﬁnancial ﬁrm is unable to function when the value of its equity falls to a sufﬁciently small fraction of its outstanding liabilities. In good times, such a ﬁrm is likely be acquired, may be able to raise new capital or may face an orderly bankruptcy. If this capital shortage occurs just when the ﬁnancial sector is already ﬁnancially constrained, then the government faces the question of whether to rescue the ﬁrm with taxpayer money as other avenues are no longer available. Such a capital shortage is damaging to the real economy as the failure of this ﬁrm will have repercussions throughout the ﬁnancial and real sectors. Consequently a ﬁrm is systemically risky if it is likely to face a large capital shortfall just when the ﬁnancial sector itself is under distress. The goal of this paper is to propose an empirical methodology to measure capital shortage in order to produce a systemic risk index for an individual ﬁrm. The capital shortfall of a ﬁrm depends on its degree of leverage and the equity loss that would result from such a crisis. While the degree of leverage can be measured, the equity loss in a crisis must be predicted. We predict it with the Marginal Expected Shortfall (MES), the expected equity loss of a ﬁrm 2

不可压缩均匀各向同性湍流直接数值 1024 POF Gotoh_Fukayama_Nakano

Velocityﬁeld statistics in homogeneous steady turbulence obtained using a high-resolution direct numerical simulationToshiyuki Gotoh a)Department of Systems Engineering,Nagoya Institute of Technology,Showa-ku,Nagoya466-8555,JapanDaigen FukayamaInformation and Mathematical Science Laboratory,Inc.,2-43-1,Ikebukuro,Toshima-ku,Tokyo171-0014,JapanTohru NakanoDepartment of Physics,Chuo University,Kasuga,Bunkyo-ku,Tokyo112-8551,JapanVelocityﬁeld statistics in the inertial to dissipation range of three-dimensional homogeneous steadyturbulentﬂow are studied using a high-resolution DNS with up to Nϭ10243grid points.The rangeof the Taylor microscale Reynolds number is between38and460.Isotropy at the small scales ofmotion is well satisﬁed from half the integral scale͑L͒down to the Kolmogorov scale͑␩͒.TheKolmogorov constant is1.64Ϯ0.04,which is close to experimentally determined values.The thirdorder moment of the longitudinal velocity difference scales as the separation distance r,and itscoefﬁcient is close to4/5.A clear inertial range is observed for moments of the velocity differenceup to the tenth order,between2␭Ϸ100␩and L/2Ϸ300␩,where␭is the Taylor microscale.Thescaling exponents are measured directly from the structure functions;the transverse scalingexponents are smaller than the longitudinal exponents when the order is greater than four.Thecrossover length of the longitudinal velocity structure function increases with the order andapproaches2␭,while that of the transverse function remains approximately constant at␭.Thecrossover length and importance of the Taylor microscale are discussed.©2002AmericanInstitute of Physics.͓DOI:10.1063/1.1448296͔I.INTRODUCTIONKolmogorov studied the statistical laws of a velocity ﬁeld for small scales of turbulent motion at high Reynolds numbers.1,2Two hypotheses were introduced in his theory ͑hereafter K41for short͒:local isotropy and homogeneity exists;and there is an inertial range in the energy spectrum of theﬂow that is independent of viscosity and large-scale properties at sufﬁciently high Reynolds numbers.The most prominent conclusion of his theory is the presence of the Kolmogorov spectrum E(k)ϭK⑀¯2/3kϪ5/3in the inertial range,where⑀¯is the average rate of energy dissipation per unit mass and K is a universal constant.Since K41,there has been a considerable amount of ef-fort made to study the turbulent velocityﬁeld statistics in the inertial range,and the energy spectrum has been a central quantity of interest.The Kolmogorov spectrum and constant have been measured inﬁeld and laboratory experiments.3–7 The exponent for the inertial range spectrum is now widely accepted asϪ5/3,with a small correction to account forﬂow intermittency.The Kolmogorov constant K is between1.5 and 2.After studying the results of many experiments, Sreenivasan stated that K is1.62Ϯ0.17.7The spectral theory of turbulence has also been used to predict the Kolmogorov constant.The value of K is1.77when the Lagrangian history direct interaction approximation͑LHDIA͒is used,8,9and 1.72when the Lagrangian renormalized approximation ͑LRA͒is used.10,11These are fully systematic theories that do not contain any ad hoc parameters.Direct numerical simulations͑DNSs͒of turbulentﬂows are now performed at higher Reynolds numbers,due to the recent dramatic increase in computational power.In the early 90’s,the resolution of DNS reached Nϭ5123grid points with a Taylor microscale Reynolds number R␭of 210ϳ240.12–20Most high-resolution DNSs have been per-formed for steady turbulence conditions to achieve high Rey-nolds numbers and obtain reliable statistics.Although results were reported with R␭greater than200,an inertial range spectrum was observed only for the lowest narrow wave number band at which forcing was applied.The Kolmogorov constant was inferred to be about1.5ϳ2in Ref.14and1.62 in Ref.16,but these results are not convincing,due to the insufﬁcient width of the scaling range,anisotropy of theﬂow ﬁeld,limited ensemble size,forcing techniques used,and numerical limitations of the simulations.Intermittency has also attracted the interest of research-ers.Since Kolmogorov’s intermittency theory͑hereafter K62͒,21many theoretical and statistical models of intermit-tency have been developed.4,22,23The scaling exponents of higher order structure functions for velocity differences in the inertial range were studied intensively.Intermittency in-creases with a decrease in the size of the scales of motion. The small-scale statistics gradually deviate from a Gaussian distribution,and the scaling exponents differ from those pre-dicted by K41.Experiments at very high Reynolds numbers have beena͒Electronic mail:gotoh@system.nitech.ac.jpPHYSICS OF FLUIDS VOLUME14,NUMBER3MARCH200210651070-6631/2002/14(3)/1065/17/$19.00©2002American Institute of Physicsperformed in the atmospheric boundary layer and in huge wind tunnels,and the measured scaling exponents were found to deviate from K41scaling.5,6,24,25,84However,there have been arguments made about the lack of small-scaleﬂow isotropy and homogeneity in these experiments,which might be affected by the large-scale shear.25,26For experiments at moderate Reynolds numbers under relatively well-controlled laboratory conditions,the width of the scaling range is usually not large enough to determine the scaling exponents precisely.Extended self-similarity͑ESS͒has been exploited to overcome this difﬁculty and applied to various turbulentﬂows in both experiments and DNSs.28–31 The idea is to measure the scaling exponents of the structure functions when they are plotted against the third order lon-gitudinal structure function,rather than to use the separation distance.The width of the scaling range is longer than that obtained with the usual method at low to moderate Reynolds numbers.The scaling exponents are anomalous,but do agree with those obtained from high Reynolds number experiments up to a certain order.24,28–31However,there is no consensus as to why the structure functions give a longer inertial range, or what is missing from theﬂow statistics as a result.Also, there is no unique way to determine the scaling exponents for the transverse and mixed velocity structure functions,be-cause those higher order structure functions can be plotted against other types of third order structure functions as well as the third order longitudinal structure function.There also have been arguments about whether the scal-ing exponents for the longitudinal and transverse structure functions at small scales are equal.25–27,32–38Many experi-ments and DNSs have reported that higher order longitudinal scaling exponents are larger than transverse ones.However, some researchers have argued that the difference is due to deviation from the assumed conditions,such as local homo-geneity,isotropy,and the independence of small scales from macroscale parameters.They have suggested that when the Reynolds number becomes large enough,the difference will vanish.36,37,39,40In many aspects of turbulence research,there have been questions posed about the extent to which the local homoge-neity and isotropy of the turbulent velocityﬁeld are attained. This will affect the small-scale statistics signiﬁcantly.Recent experimental studies have shown that local isotropy is par-tially satisﬁed for lower order moments.25,26,37However,it is not sufﬁcient to examine only the conditions assumed in the above studies,and only a limited knowledge of the trueﬂow conditions is available so far.26,37,38A DNS with a sufﬁciently large grid size provides abetter opportunity to examine the points raised above.It hasthe advantage that any physical quantity can be measureddirectly without deforming theﬂowﬁeld.In the presentstudy,a series of large scale DNSs have been performed at ahigh resolution of up to Nϭ10243and R␭ϭ460.41–45The inertial range of the turbulenceﬁeld has a considerablelength,and useful velocity statistics can be extracted such asthe Kolmogorov constant,the energy spectrum,velocitystructure functions up to the tenth order,their scaling expo-nents,and probability density functions for velocity differ-ences.To the authors’best knowledge,these are theﬁrstDNS data in the inertial range;the data provide new insightinto the inertial and dissipation ranges.The main purposes of the present paper are to describethe statistics of the velocityﬁeld in an incompressible steadyturbulentﬂow obtained from the DNS,and to reexaminecurrent knowledge of turbulence,developed since K41.Thepaper is organized as follows.The numerical aspects of thepresent DNS are described in Sec.II,and the energy spec-trum is examined in Sec.III.The variation of single pointquantities and probability density functions͑PDFs͒with theReynolds number is discussed in Sec.IV.The isotropy of thesecond and third order moments of the velocity difference isexamined in Sec.V,and the energy budget is examined interms of the Ka´rma´n–Howarth–Kolmogorov equation inSec.VI.The structure functions and scaling exponents arediscussed in Sec.VII.Section VIII presents an analysis ofthe crossover lengths of the structure functions.Finally,asummary and conclusions are provided in Sec.IX.II.NUMERICAL SIMULATIONThe Navier–Stokes equations are integrated in Fourierspace for unit density:ͩץץtϩ␯k2ͪuϭP͑k͒•F͓uÃ␻͔kϩf,͑1͒͗f͑k,t͒f͑Ϫk,s͒͘ϭP͑k͒F͑k͒4␲k2␦͑tϪs͒,͑2͒where␻is the vorticity vector,P͑k͒is the projection opera-tor,F denotes a Fourier transform,and f is a solenoidal Gaussian random force that is white in time.The spectrum of the random force F(k)is constant over the low wave number band and zero otherwise;the force is normalized asTABLE I.DNS parameters and statistical quantities of the runs.T eddya v is the period used for the time average.R␭N k max␯c f Forcing range T eddya v E⑀¯L␭␩(ϫ10Ϫ2)K 38128360 1.50ϫ10Ϫ2 1.30ͱ3рkрͱ1222.6 1.99 1.190.8910.501 4.10¯5425631217.00ϫ10Ϫ30.70ͱ3рkрͱ1214.9 1.390.6270.8290.393 2.72¯702563121 4.00ϫ10Ϫ30.50ͱ3рkрͱ1249.7 1.160.4570.7850.318 1.93¯1255123241 1.35ϫ10Ϫ30.50ͱ3рkрͱ12 5.52 1.250.4920.7440.1850.841¯2845123241 6.00ϫ10Ϫ40.501рkрͱ6 3.03 1.960.530 1.2460.1490.449 1.64 38110243483 2.80ϫ10Ϫ40.511рkрͱ6 4.21 1.740.499 1.1390.09890.258 1.63 46010243483 2.00ϫ10Ϫ40.511рkрͱ6 2.14 1.790.506 1.1500.08410.199 1.64 1066Phys.Fluids,Vol.14,No.3,March2002Gotoh,Fukayama,and Nakano͵ϱF ͑k ͒dk ϭ⑀¯in ,͑3͒where ⑀¯in is the average rate of the energy input per unit mass.A pseudo-spectral code was used to compute the con-volution sums,and the aliasing error was effectively re-moved.The time integration was performed using the fourth order Runge–Kutta–Gill method.Physical quantities of turbulent ﬂow include the total energyE ͑t ͒ϭ12͗u 2͘ϭ32u ¯2ϭ͵ϱE ͑k ͒dk ,͑4͒the average energy dissipation per unit mass⑀¯ϭ2␯͵ϱk 2E ͑k ͒dk ,͑5͒the integral scaleL ϭͩ3␲4͵ϱk Ϫ1E ͑k ͒dkͪͲE ,͑6͒the Taylor microscale␭ϭͩ5EͲ͵ϱk 2E ͑k ͒dkͪ1/2,͑7͒the Taylor microscale Reynolds numberR ␭ϭu ¯␭␯,͑8͒and the Kolmogorov scale␩ϭͩ␯3⑀¯ͪ1/4.͑9͒The range of the Taylor microscale Reynolds number was 38to 460.The characteristic parameters of the DNS are listed in Table I.43Most of these are identical to Gotoh and Fukayama,43but the averaging time for R ␭ϭ381was ex-tended to 4.21large eddy turnover times.A statistically steady state was conﬁrmed by observing the time evolutionof the total energy,the total enstrophy,and the skewness of the longitudinal velocity derivative.The statistical averages were computed as time averages over tens of large eddy turnover times for the lower Reynolds number ﬂows,and over a few large eddy turnover times for the higher Reynolds number ﬂows.The resolution condition k max ␩Ͼ1was satis-ﬁed for most runs,except for R ␭ϭ460in which k max ␩was slightly less than unity (k max ␩ϭ0.96).This does not ad-versely affect the results in the inertial range.The computational time required for runs at a N ϭ10243resolution varied,depending on the statistical data that was gathered.Typically,60h was required for one large eddy turnover time.The total time of the computations was more than 500h for the longest run (R ␭ϭ381).Data col-lected during the transition period to steady state ͑about six large eddy turnover times ͒were discarded.The relatively long time required to attain steady state was due to the low wave number band forcing.This imposes a severe computa-tional putations with R ␭р284were per-formed on a Fujitsu VPP700E parallel vector machine with 16processors at RIKEN.Simulations of higher R ␭were per-formed on a Fujitsu VPP5000/56with 32processors at the Nagoya University Computation Center.III.ENERGY SPECTRUMFigure 1shows the three-dimensional energy spectrum calculated for each run.All of the curves are scaled to the Kolmogorov units and multiplied by k 5/3.As the Reynolds number increases,the curves extend toward lower wave numbers.The curves of ﬂows with Reynolds numbers larger than R ␭ϭ284contain a ﬁnite plateau,which indicates that E (k )ϰk Ϫ5/3.There is a bump when 0.04рk ␩р0.3at the high end of the inertial range,which is consistent with pre-vious experimental and numerical observations.6,16The nor-malized energy transfer ﬂux,deﬁned by1⑀¯⌸͑k ͒ϭ1⑀¯͵kϱT ͑k Ј͒dk Ј͑10͒is shown in Fig.2,where T (k )is a nonlinear energy transfer function in the energy spectrum equation.4,22Between0.007рk ␩р0.04,⌸(k )/⑀¯is approximately constantand FIG.1.Scaled energy spectra,⑀¯Ϫ1/4␯Ϫ5/4(k ␩)5/3E (k ).The inertial range is 0.007рk ␩р0.04and K ϭ1.64Ϯ0.04.A horizontal line indicates K ϭ1.64.FIG.2.Normalized energy transfer ﬂux,⌸(k )/⑀¯for R ␭ϭ381and 460.1067Phys.Fluids,Vol.14,No.3,March 2002Velocity ﬁeld statistics in homogeneousclose to unity;thus the ﬂow is in an equilibrium state over the inertial range of the energy spectrum,corresponding to the plateaus in Fig.1.The Kolmogorov constant given in Table I is determined using a least square ﬁt between 0.007рk ␩р0.04on the R ␭Ͼ284curves.In Ref.43,the Kolmog-orov constant was reported as K ϭ1.65Ϯ0.05.However,the averaging time has since been extended for the R ␭ϭ381run.The R ␭ϭ478run differs slightly from statisticalequilibrium,since ⌸(k )/⑀¯is not exactly one;for this reason,the R ␭ϭ478data were not used for this analysis.The Kol-mogorov constant,computed using the data only from the R ␭ϭ381and 460runs,isK ϭ1.64Ϯ0.04,͑11͒which is in good agreement with experimental values and recent DNS data.7,16There are many DNSs reporting the Kolmogorov constant higher than the value 1.64.However,the length of the inertial range in those DNSs is not long enough to clearly observe the k Ϫ5/3range,and the top of the bump of the compensated energy spectrum k 5/3E (k )is un-derstood as the inertial range,so that the Kolmogorov con-stant is read as about 2as seen in Fig.1.16The Kolmogorov constant 1.64is also close to the value obtained using the LHDIA ͑1.77͒,8,9the LRA ͑1.72͒.10,11These spectral theories of turbulence are consistent with Lagrangian dynamics,arederived systematically,and contain no ad hoc parameters.Figure 3shows the one-dimensional energy spectrum ob-tained from the present DNS with R ␭ϭ460,from experi-ments,and from the LRA.The agreement between the curves is satisfactory.Therefore we conclude that the present DNS has successfully calculated a homogeneous turbulent ﬂow ﬁeld in the inertial range of the energy spectrum.IV.ONE-POINT STATISTICS A.MomentsSome one-point moments of the velocity ﬁeld areS 3͑u ͒ϵ͗u 3͗͘u 2͘3/2,S 3͑u x ͒ϵ͗u x 3͗͘x 2͘3/2,͑12͒K 4͑u ͒ϵ͗u 4͗͘u 2͘2,K 4͑u x ͒ϵ͗u x 4͗͘u x 2͘2,K 4͑u y ͒ϵ͗u y 4͗͘u y 2͘2,͑13͒where u is the velocity component in the x direction.The variation of these moments with the Reynolds number is shown in Fig.4and listed in Table II.The general behavior of the curves is consistent with previous DNS and experi-mental data.13,14,18,19,26,46,47There are small effects of rela-tively low resolution on S 3and K 4for the velocity deriva-tives for R ␭ϭ381and 460data.The skewness factor of the velocity u is very small for runs with the R ␭р125,and isofparison of one-dimensional energy spectra.Symbols:experi-ments,solid line:present DNS (R ␭ϭ460),dashed line:statistical theory ͑LRA and MLRA ͒.FIG.4.Variation of the moments of the velocity and velocity gradient withthe Reynolds number.Line:present DNS,circle:K 4(u y )͑Jime´nez et al.,Ref.13͒,solid square:K 4(u )͑Jime ´nez et al.,Ref.13͒,square:K 4(u x )͑Wang et al.,Ref.14͒,plus:K 4(u x )͑Vedula and Yeung,Ref.18͒,star:ϪS 3(u x )͑Wang et al.,Ref.14͒.TABLE II.Moments of the velocity and velocity derivatives.R ␭S 3(u )K 4(u )S 3(ץu /ץx )K 4(ץu /ץx )K 4(ץu /ץy )380.0227 2.89Ϫ0.520 4.14 5.16540.00563 2.86Ϫ0.517 4.47 6.00700.00473 2.93Ϫ0.519 4.81 6.621250.0820 2.94Ϫ0.529 5.658.192840.0231 2.77Ϫ0.531 6.6310.1381Ϫ0.246 2.98Ϫ0.5747.9012.2460Ϫ0.1682.89Ϫ0.5457.9111.71068Phys.Fluids,Vol.14,No.3,March 2002Gotoh,Fukayama,and Nakanothe order of 0.2for runs with the R ␭у284.The relatively large values of the velocity skewness are caused by the shorter averaging time used compared to the low Reynolds number runs.Since most of the energy resides in the lowest wave number band,there are persistent large ﬂuctuations of the large scales of motion over longer time period.The longer time average or the forcing at larger wave numbers would yield smaller velocity skewness.The ﬂatness factor of the velocity ﬁeld is close to three,which is the Gaussian value.The skewness factor of the longitudinal velocity deriva-tives is very insensitive to the Reynolds number,S 3͑u x ͒ϰR ␭0.0370,͑14͒where the exponent is determined by a least square ﬁt.Theaverage value is Ϫ0.53,which is consistent with experimen-tal observations over the range of Reynolds numbers studied in the present work.However,the exponent is smaller than indicated by the experimental data.26,46The ﬂatness factors for the longitudinal and transverse velocity derivatives in-crease with the Reynolds number asK 4͑u x ͒ϰR ␭0.266,K 4͑u y ͒ϰR ␭0.335.͑15͒The exponent of K 4(u y )is larger than that of K 4(u x );thus,the PDF for the transverse velocity derivative has longer tails than those of the longitudinal velocity derivative.From ex-perimental observations,Shen and Warhaft reported thatK 4(u x )ϰR ␭0.37and K 4(u y )ϰR ␭0.25.26Since there is scatter in the experimetal data,the exponents in Eq.͑15͒by the present DNS are not inconsistent with the experimental data.Van Atta and Antonia studied the Reynolds number dependence of S 3(u x )and K 4(u x ),46and found thatS 3͑u x ͒ϰR ␭0.12,K 4͑u x ͒ϰR ␭0.32for ␮ϭ0.2,͑16͒S 3͑u x ͒ϰR ␭0.15,K 4͑u x ͒ϰR ␭0.41for ␮ϭ0.25,͑17͒where ␮is the exponent deﬁned by ͗⑀r 2͘ϰr Ϫ␮for the locally averaged energy dissipation rate.4,21Generally,the Reynolds number dependency of S 3and K 4in our DNSs is weaker than observed in the experiments,irrespective of the type of forcing used.We believe this is because the range of Rey-nolds numbers in DNS is smaller than experimental ﬂows,and there remain small-scale anisotropy effects in the experi-ments.B.Probability density functionsThe probability density function conveys information about single-point velocity statistics.It has been one of the central issues of turbulence research in the last decade.Single-point PDFs for the velocity and its derivatives are shown in Figs.5–7.A longer time period was necessary for the time average to obtain well-converged PDF for the ve-locity Q (u ).The distribution Q (u )is close to Gaussian,and its tail extends to very low values of the order of 10Ϫ10.Such values have not been reported in the literature.The Q (u )curve for R ␭ϭ381is skewed negatively,but this is attributed to the insufﬁcient time-averaging period ͑four large eddy turnover times ͒that was used.The overall trend is that Q (u )decays faster than a Gaussian distribution at large ampli-tudes.This behavior was also observed in one-dimensional decaying and forced Burgers turbulence.48,49Jime´nez has shown that the PDF Q (u )is slightly sub-Gaussian as the energy spectrum decays faster than k Ϫ1.50FIG.5.Variation of velocity PDF with the Reynoldsnumber.FIG.6.Variation of the longitudinal velocity derivative PDF with the Rey-noldsnumber.FIG.7.Variation of the transverse velocity derivative PDF with the Rey-nolds number.1069Phys.Fluids,Vol.14,No.3,March 2002Velocity ﬁeld statistics in homogeneousThis is consistent with the present DNS results.Studies of the Q (u )tail predict that Q (w )ϰexp(Ϫc ͉w ͉3)when the forc-ing has a short correlation time.51,52Here,w ϭu /͗u 2͘1/2is the normalized velocity amplitude and c is a nondimensional constant.The asymptotic form of Q (u )was examined by plotting ln ͓Ϫln(Q (w )͔against ln ͉w ͉;however,the Q (w )tails were too short to determine the true asymptotic form.The PDF for the longitudinal velocity derivative is slightly skewed,as expected from the ﬁnite negative value of the skewness factor.The tail becomes longer as the Reynolds number increases.Figure 7shows that the PDF of the trans-verse derivative is symmetric and has a longer tail than the longitudinal derivative.There are many theories for the PDF of the velocity derivative.The asymptotic tail of Q (ץu /ץy )is presented in Fig.8,in which both the positive and negative sides are plotted by assuming that the PDF is symmetric.The tails gradually become longer as the Reynolds number increases;therefore,Q (s )is Reynolds-number dependent,and cannot be represented in a single stretched exponential form as Q (s )ϰexp(Ϫb ͉s ͉h ),where s is the normalized amplitude of ץu /ץy and b is a nondimensional constant that is a function of the Reynolds number.53V.ISOTROPYThe hypothesis of isotropy of the ﬂow ﬁeld is one of the key components of K41.There are various methods to ex-amine the degree of isotropy.One measure of isotropy can be obtained from the relations between the second and third order longitudinal and transverse velocity structure func-tions.These areD LL ϵ͗͑␦u r ͒2͘,D TT ϵ͗͑␦v r ͒2͘,͑18͒D LLL ϵ͗͑␦u r ͒3͘,D LTT ϵ͗␦u r ͑␦v r ͒2͘,͑19͒where␦u r ϵ͑u ͑x ϩr ͒Ϫu ͑x ͒͒•r /r ,͑20͒␦v r ϵ͑u ͑x ϩr ͒Ϫu ͑x ͒͒•͑I Ϫrr /r 2͒•e Ќ,͑21͒and e Ќis the unit vector perpendiculer to r ,and I is the unit tensor.Then the isotropy and incompressibility relations areD TT ͑r ͒ϭD LL ͑r ͒ϩr 2dD LL ͑r ͒dr ,͑22͒D LTT ͑r ͒ϭ16ddrrD LLL ͑r ͒.͑23͒In DNS,the solenoidal property of the Fourier amplitude velocity vector u ͑k ͒is always satisﬁed to the level of nu-merical error,which is smaller than 10Ϫ15.Thus,the accu-racy of the above relations depends solely on the deviation from isotropy.The two sides of Eqs.͑22͒and ͑23͒are com-pared for R ␭ϭ125,381,and 460in Figs.9and 10.The curves in the ﬁgures are divided by r 2/3and r ,respectively,and the vertical axes of the plots are linear.The thick lines represent the left hand sides of Eqs.͑22͒and ͑23͒,and the thin lines correspond to the right-hand sides.The isotropy of the second and third order moments is excellent for scales less than L /2.The difference at larger separations is caused by the anisotropy due to the small number of energy-containing Fourier modes.The curves for R ␭ϭ381and460FIG.8.Variation of the asymptotic tail of the transverse velocity derivative PDF with the Reynolds number.Both positive and negative sides are plot-ted.The rightmost curve corresponds to R ␭ϭ460.FIG.9.Isotropy relation at the second order.Thin line:D TT (r )r Ϫ2/3,thick line:(D LL (r )ϩ(r /2)(dD LL (r )/dr ))r Ϫ2/3.L /␩and ␭/␩are shown for R ␭ϭ460.FIG.10.Isotropy relation at the third order.Thin line:D LTT (r )r Ϫ1,thick line:((1/6)(d /dr )rD LLL (r ))r Ϫ1.L /␩and ␭/␩are shown for R ␭ϭ460.1070Phys.Fluids,Vol.14,No.3,March 2002Gotoh,Fukayama,and Nakanoin Fig.9are not horizontal,suggesting that the second order structure function does not scale as r 2/3.The scaling expo-nents will be examined later in this paper.The isotropic re-lations,such as D 1122ϭD 1133and D 2222ϭ3D 2233ϭD 3333,and Hill’s higher order relations were not computed.54VI.KA´RMA ´N–HOWARTH–KOLMOGOROV EQUATION The energy budget for various scales is described by the Ka´rma ´n–Howarth–Kolmogorov ͑KHK ͒equation,45⑀¯r ϭϪD LLL ϩ6␯ץD LL ץrϩZ ͑24͒for steady turbulence,4,55,56where Z (r )denotes contributions due to the external force given by Z ͑r ,t ͒ϭ͵Ϫϱt͗␦f ͑r ,t ͒•␦f ͑r ,s ͒͘dsϭ12r͵0ϱͩ115ϩsin kr ͑kr ͒3ϩ3cos kr ͑kr ͒4Ϫ3sin kr ͑kr ͒5ͪF ͑k ͒dk .͑25͒Since the external force spectrum F (k )is localized in arange of low wave numbers,the asymptotic form of Z (r )for small separations is given asZ ͑r ͒ϭ235⑀¯in k f 2r 3,k f 2ϵ͐0ϱk 2F ͑k ͒dk͐0ϱF ͑k ͒dk.͑26͒A generalized Ka´rma ´n–Howarth–Kolmogorov equation has also been derived:57–6343⑀¯r ϭϪ͑D LLL ϩ2D LTT ͒ϩ2␯ץץr͑D LL ϩ2D TT ͒ϩW ,͑27͒where W ͑r ͒ϭ4r͵0ϱͩ13ϩcos kr ͑kr ͒2Ϫsin kr͑kr ͒3ͪF ͑k ͒dk ,Ϸ215⑀¯in r 3k f 2for ͉k f r ͉Ӷ1.͑28͒Equation ͑24͒is recovered by substituting Eqs.͑22͒and ͑23͒into Eq.͑27͒.Figure 11shows the results obtained when each term of Eq.͑24͒is divided by ⑀¯r for R ␭ϭ460.Curves in which r /␩is larger than r /␩ϭ1200are not shown,because the sign of D LLL changes.A thin horizontal line indicates the Kolmog-orov value 4/5.When the separation distance decreases,the effect of the large scale forcing used in the present DNS decreases quickly,while the viscous term grows gradually.The third order longitudinal structure function D LLL quickly rises to the Kolmogorov value,remains there over the iner-tial range ͑between r /␩Ϸ50and 300͒,and then decreases.In the inertial range,the force term decreases as r 3according to Eq.͑26͒,while the viscous term increases as r ␨2Ϫ1(␨2Ͻ1)when r decreases.͓Since each term in the ﬁgure is divided by (⑀¯r ),the slope of each curve is 2and ␨2Ϫ2,respectively.͔The sum of the three terms in the right hand side of Eq.͑24͒divided by ⑀¯r is close to 4/5,the Kolmogorov value.The deviation of the sum from the 4/5law at the smallest scales is due to the slightly lower resolution of the data at these scales ͑k max ␩is close to one ͒.At larger scales greater than r /␩ϭ700,the deviation is caused by the ﬁniteness oftheFIG.11.Terms in the Ka´rma ´n–Howarth–Kolmogorov equation when R ␭ϭ460.Thin solid line:4/5.FIG.12.Kolmogorov’s 4/5law.L /␩and ␭/␩are shown for R ␭ϭ460.The maximum values of the curves are 0.665,0.771,0.781,and 0.757for R ␭ϭ125,284,381,and 460,respectively.FIG.13.Terms in the generalized Ka´rma ´n–Howarth–Kolmogorov equation for R ␭ϭ460.Thin solid line:4/3.1071Phys.Fluids,Vol.14,No.3,March 2002Velocity ﬁeld statistics in homogeneousensemble,which indicates the persistent anisotropy of the larger scales.The above ﬁndings are consistent with the cur-rent knowledge of turbulence developed since Kolmogorov,although conﬁrmation of some aspects of turbulence using actual data is new from both a numerical and experimental point of view.56,59–65It is interesting and important to observe when the Kol-mogorov 4/5law is satisﬁed as the Reynolds numberincreases.6,66–69Figure 12shows curves of ϪD LLL (r )/(⑀¯r )for various Reynolds numbers.In this ﬁgure,the 4/5law applies when the curves are horizontal.The portion of the curves in which r /␩Ͼ1200is not shown.Although there is a small but ﬁnite horizontal range when R ␭Ͼ284,the level of the plateau is still less than the Kolmogorov value.The maximum values of the curves are 0.665,0.771,0.781,and 0.757for R ␭ϭ125,284,381,and 460,respectively.The value 0.781for R ␭ϭ381is 2.5%less than 0.8.An asymptotic state is approached slowly,which is consistent with recent studies.However,the asymptote is approached faster than predicted by the theoretical estimate.66,69Theslow approach is due to the fact that D LLL (r )is the third order structure function and most positive contributions are canceled by negative ones.Thus only the slight asymmetry of the ␦u r PDF contributes to D LLL .The level of the plateau of the R ␭ϭ460curve is slightly less than the others.A higher value would be expected if the time average period used for the R ␭ϭ460run were longer.The generalized Ka´rma ´n–Howarth–Kolmogorov equa-tion Eq.͑27͒is also examined in a similar fashion.Figure 13shows each term of the equation divided by ⑀¯r ;a horizontal line indicates the 4/3law.The agreement between the present data and theory is satisfactory.The third order moment slowly approaches the Kolmogorov value 4/3,as shown in Fig.14.The maximum values of the curves of the 4/3law are 0.564,1.313,1.297,and 1.259for R ␭ϭ125,284,381,and 460,respectively.VII.STRUCTURE FUNCTIONS AND SCALING EXPONENTSThe velocity structure functions are deﬁned asS p L ͑r ͒ϭ͉͗␦u r ͉p͘,S p T ͑r ͒ϭ͉͗␦v r ͉p͘,FIG.14.Kolmogorov’s 4/3law.L /␩and ␭/␩are shown for R ␭ϭ460.The maximum values of the curves for the 4/3law are 0.564,1.313,1.297,and 1.259for R ␭ϭ125,284,381,and 460,respectively.FIG.15.Variation of the ␦u r PDF with r for R ␭ϭ381.From the outermostcurve,r n /␩ϭ2n Ϫ1dx /␩ϭ2.38ϫ2n Ϫ1,n ϭ1,...,10,where dx ϭ2␲/1024.The inertial range corresponds to n ϭ6,7,8.Dotted line:Gaussian.FIG.16.Variation of PDF for ␦v r with r at R ␭ϭ381.The classiﬁcation of curves is the same as in Fig.17.FIG.17.Convergence of the tenth order accumulated moments C 10(␦u r )at R ␭ϭ381for various separations r n /␩ϭ2.38ϫ2n Ϫ1,n ϭ1,...,10.Curves are for n ϭ1,...10from the uppermost,and the inertial range corresponds to n ϭ6,7,8.1072Phys.Fluids,Vol.14,No.3,March 2002Gotoh,Fukayama,and Nakano。

T.W. ANDERSON (1971). The Statistical Analysis of Time Series. Series in Probability and Ma

425 BibliographyH.A KAIKE(1974).Markovian representation of stochastic processes and its application to the analysis of autoregressive moving average processes.Annals Institute Statistical Mathematics,vol.26,pp.363-387. B.D.O.A NDERSON and J.B.M OORE(1979).Optimal rmation and System Sciences Series, Prentice Hall,Englewood Cliffs,NJ.T.W.A NDERSON(1971).The Statistical Analysis of Time Series.Series in Probability and Mathematical Statistics,Wiley,New York.R.A NDRE-O BRECHT(1988).A new statistical approach for the automatic segmentation of continuous speech signals.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-36,no1,pp.29-40.R.A NDRE-O BRECHT(1990).Reconnaissance automatique de parole`a partir de segments acoustiques et de mod`e les de Markov cach´e s.Proc.Journ´e es Etude de la Parole,Montr´e al,May1990(in French).R.A NDRE-O BRECHT and H.Y.S U(1988).Three acoustic labellings for phoneme based continuous speech recognition.Proc.Speech’88,Edinburgh,UK,pp.943-950.U.A PPEL and A.VON B RANDT(1983).Adaptive sequential segmentation of piecewise stationary time rmation Sciences,vol.29,no1,pp.27-56.L.A.A ROIAN and H.L EVENE(1950).The effectiveness of quality control procedures.Jal American Statis-tical Association,vol.45,pp.520-529.K.J.A STR¨OM and B.W ITTENMARK(1984).Computer Controlled Systems:Theory and rma-tion and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.M.B AGSHAW and R.A.J OHNSON(1975a).The effect of serial correlation on the performance of CUSUM tests-Part II.Technometrics,vol.17,no1,pp.73-80.M.B AGSHAW and R.A.J OHNSON(1975b).The inﬂuence of reference values and estimated variance on the ARL of CUSUM tests.Jal Royal Statistical Society,vol.37(B),no3,pp.413-420.M.B AGSHAW and R.A.J OHNSON(1977).Sequential procedures for detecting parameter changes in a time-series model.Jal American Statistical Association,vol.72,no359,pp.593-597.R.K.B ANSAL and P.P APANTONI-K AZAKOS(1986).An algorithm for detecting a change in a stochastic process.IEEE rmation Theory,vol.IT-32,no2,pp.227-235.G.A.B ARNARD(1959).Control charts and stochastic processes.Jal Royal Statistical Society,vol.B.21, pp.239-271.A.E.B ASHARINOV andB.S.F LEISHMAN(1962).Methods of the statistical sequential analysis and their radiotechnical applications.Sovetskoe Radio,Moscow(in Russian).M.B ASSEVILLE(1978).D´e viations par rapport au maximum:formules d’arrˆe t et martingales associ´e es. Compte-rendus du S´e minaire de Probabilit´e s,Universit´e de Rennes I.M.B ASSEVILLE(1981).Edge detection using sequential methods for change in level-Part II:Sequential detection of change in mean.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-29,no1,pp.32-50.426B IBLIOGRAPHY M.B ASSEVILLE(1982).A survey of statistical failure detection techniques.In Contribution`a la D´e tectionS´e quentielle de Ruptures de Mod`e les Statistiques,Th`e se d’Etat,Universit´e de Rennes I,France(in English). M.B ASSEVILLE(1986).The two-models approach for the on-line detection of changes in AR processes. In Detection of Abrupt Changes in Signals and Dynamical Systems(M.Basseville,A.Benveniste,eds.). Lecture Notes in Control and Information Sciences,LNCIS77,Springer,New York,pp.169-215.M.B ASSEVILLE(1988).Detecting changes in signals and systems-A survey.Automatica,vol.24,pp.309-326.M.B ASSEVILLE(1989).Distance measures for signal processing and pattern recognition.Signal Process-ing,vol.18,pp.349-369.M.B ASSEVILLE and A.B ENVENISTE(1983a).Design and comparative study of some sequential jump detection algorithms for digital signals.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-31, no3,pp.521-535.M.B ASSEVILLE and A.B ENVENISTE(1983b).Sequential detection of abrupt changes in spectral charac-teristics of digital signals.IEEE rmation Theory,vol.IT-29,no5,pp.709-724.M.B ASSEVILLE and A.B ENVENISTE,eds.(1986).Detection of Abrupt Changes in Signals and Dynamical Systems.Lecture Notes in Control and Information Sciences,LNCIS77,Springer,New York.M.B ASSEVILLE and I.N IKIFOROV(1991).A uniﬁed framework for statistical change detection.Proc.30th IEEE Conference on Decision and Control,Brighton,UK.M.B ASSEVILLE,B.E SPIAU and J.G ASNIER(1981).Edge detection using sequential methods for change in level-Part I:A sequential edge detection algorithm.IEEE Trans.Acoustics,Speech,Signal Processing, vol.ASSP-29,no1,pp.24-31.M.B ASSEVILLE, A.B ENVENISTE and G.M OUSTAKIDES(1986).Detection and diagnosis of abrupt changes in modal characteristics of nonstationary digital signals.IEEE rmation Theory,vol.IT-32,no3,pp.412-417.M.B ASSEVILLE,A.B ENVENISTE,G.M OUSTAKIDES and A.R OUG´E E(1987a).Detection and diagnosis of changes in the eigenstructure of nonstationary multivariable systems.Automatica,vol.23,no3,pp.479-489. M.B ASSEVILLE,A.B ENVENISTE,G.M OUSTAKIDES and A.R OUG´E E(1987b).Optimal sensor location for detecting changes in dynamical behavior.IEEE Trans.Automatic Control,vol.AC-32,no12,pp.1067-1075.M.B ASSEVILLE,A.B ENVENISTE,B.G ACH-D EVAUCHELLE,M.G OURSAT,D.B ONNECASE,P.D OREY, M.P REVOSTO and M.O LAGNON(1993).Damage monitoring in vibration mechanics:issues in diagnos-tics and predictive maintenance.Mechanical Systems and Signal Processing,vol.7,no5,pp.401-423.R.V.B EARD(1971).Failure Accommodation in Linear Systems through Self-reorganization.Ph.D.Thesis, Dept.Aeronautics and Astronautics,MIT,Cambridge,MA.A.B ENVENISTE and J.J.F UCHS(1985).Single sample modal identiﬁcation of a nonstationary stochastic process.IEEE Trans.Automatic Control,vol.AC-30,no1,pp.66-74.A.B ENVENISTE,M.B ASSEVILLE and G.M OUSTAKIDES(1987).The asymptotic local approach to change detection and model validation.IEEE Trans.Automatic Control,vol.AC-32,no7,pp.583-592.A.B ENVENISTE,M.M ETIVIER and P.P RIOURET(1990).Adaptive Algorithms and Stochastic Approxima-tions.Series on Applications of Mathematics,(A.V.Balakrishnan,I.Karatzas,M.Yor,eds.).Springer,New York.A.B ENVENISTE,M.B ASSEVILLE,L.E L G HAOUI,R.N IKOUKHAH and A.S.W ILLSKY(1992).An optimum robust approach to statistical failure detection and identiﬁcation.IFAC World Conference,Sydney, July1993.B IBLIOGRAPHY427 R.H.B ERK(1973).Some asymptotic aspects of sequential analysis.Annals Statistics,vol.1,no6,pp.1126-1138.R.H.B ERK(1975).Locally most powerful sequential test.Annals Statistics,vol.3,no2,pp.373-381.P.B ILLINGSLEY(1968).Convergence of Probability Measures.Wiley,New York.A.F.B ISSELL(1969).Cusum techniques for quality control.Applied Statistics,vol.18,pp.1-30.M.E.B IVAIKOV(1991).Control of the sample size for recursive estimation of parameters subject to abrupt changes.Automation and Remote Control,no9,pp.96-103.R.E.B LAHUT(1987).Principles and Practice of Information Theory.Addison-Wesley,Reading,MA.I.F.B LAKE and W.C.L INDSEY(1973).Level-crossing problems for random processes.IEEE r-mation Theory,vol.IT-19,no3,pp.295-315.G.B ODENSTEIN and H.M.P RAETORIUS(1977).Feature extraction from the encephalogram by adaptive segmentation.Proc.IEEE,vol.65,pp.642-652.T.B OHLIN(1977).Analysis of EEG signals with changing spectra using a short word Kalman estimator. Mathematical Biosciences,vol.35,pp.221-259.W.B¨OHM and P.H ACKL(1990).Improved bounds for the average run length of control charts based on ﬁnite weighted sums.Annals Statistics,vol.18,no4,pp.1895-1899.T.B OJDECKI and J.H OSZA(1984).On a generalized disorder problem.Stochastic Processes and their Applications,vol.18,pp.349-359.L.I.B ORODKIN and V.V.M OTTL’(1976).Algorithm forﬁnding the jump times of random process equation parameters.Automation and Remote Control,vol.37,no6,Part1,pp.23-32.A.A.B OROVKOV(1984).Theory of Mathematical Statistics-Estimation and Hypotheses Testing,Naouka, Moscow(in Russian).Translated in French under the title Statistique Math´e matique-Estimation et Tests d’Hypoth`e ses,Mir,Paris,1987.G.E.P.B OX and G.M.J ENKINS(1970).Time Series Analysis,Forecasting and Control.Series in Time Series Analysis,Holden-Day,San Francisco.A.VON B RANDT(1983).Detecting and estimating parameters jumps using ladder algorithms and likelihood ratio test.Proc.ICASSP,Boston,MA,pp.1017-1020.A.VON B RANDT(1984).Modellierung von Signalen mit Sprunghaft Ver¨a nderlichem Leistungsspektrum durch Adaptive Segmentierung.Doctor-Engineer Dissertation,M¨u nchen,RFA(in German).S.B RAUN,ed.(1986).Mechanical Signature Analysis-Theory and Applications.Academic Press,London. L.B REIMAN(1968).Probability.Series in Statistics,Addison-Wesley,Reading,MA.G.S.B RITOV and L.A.M IRONOVSKI(1972).Diagnostics of linear systems of automatic regulation.Tekh. Kibernetics,vol.1,pp.76-83.B.E.B RODSKIY and B.S.D ARKHOVSKIY(1992).Nonparametric Methods in Change-point Problems. Kluwer Academic,Boston.L.D.B ROEMELING(1982).Jal Econometrics,vol.19,Special issue on structural change in Econometrics. L.D.B ROEMELING and H.T SURUMI(1987).Econometrics and Structural Change.Dekker,New York. D.B ROOK and D.A.E VANS(1972).An approach to the probability distribution of Cusum run length. Biometrika,vol.59,pp.539-550.J.B RUNET,D.J AUME,M.L ABARR`E RE,A.R AULT and M.V ERG´E(1990).D´e tection et Diagnostic de Pannes.Trait´e des Nouvelles Technologies,S´e rie Diagnostic et Maintenance,Herm`e s,Paris(in French).428B IBLIOGRAPHY S.P.B RUZZONE and M.K AVEH(1984).Information tradeoffs in using the sample autocorrelation function in ARMA parameter estimation.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-32,no4, pp.701-715.A.K.C AGLAYAN(1980).Necessary and sufﬁcient conditions for detectability of jumps in linear systems. IEEE Trans.Automatic Control,vol.AC-25,no4,pp.833-834.A.K.C AGLAYAN and R.E.L ANCRAFT(1983).Reinitialization issues in fault tolerant systems.Proc.Amer-ican Control Conf.,pp.952-955.A.K.C AGLAYAN,S.M.A LLEN and K.W EHMULLER(1988).Evaluation of a second generation reconﬁgu-ration strategy for aircraftﬂight control systems subjected to actuator failure/surface damage.Proc.National Aerospace and Electronic Conference,Dayton,OH.P.E.C AINES(1988).Linear Stochastic Systems.Series in Probability and Mathematical Statistics,Wiley, New York.M.J.C HEN and J.P.N ORTON(1987).Estimation techniques for tracking rapid parameter changes.Intern. Jal Control,vol.45,no4,pp.1387-1398.W.K.C HIU(1974).The economic design of cusum charts for controlling normal mean.Applied Statistics, vol.23,no3,pp.420-433.E.Y.C HOW(1980).A Failure Detection System Design Methodology.Ph.D.Thesis,M.I.T.,L.I.D.S.,Cam-bridge,MA.E.Y.C HOW and A.S.W ILLSKY(1984).Analytical redundancy and the design of robust failure detection systems.IEEE Trans.Automatic Control,vol.AC-29,no3,pp.689-691.Y.S.C HOW,H.R OBBINS and D.S IEGMUND(1971).Great Expectations:The Theory of Optimal Stop-ping.Houghton-Mifﬂin,Boston.R.N.C LARK,D.C.F OSTH and V.M.W ALTON(1975).Detection of instrument malfunctions in control systems.IEEE Trans.Aerospace Electronic Systems,vol.AES-11,pp.465-473.A.C OHEN(1987).Biomedical Signal Processing-vol.1:Time and Frequency Domain Analysis;vol.2: Compression and Automatic Recognition.CRC Press,Boca Raton,FL.J.C ORGE and F.P UECH(1986).Analyse du rythme cardiaque foetal par des m´e thodes de d´e tection de ruptures.Proc.7th INRIA Int.Conf.Analysis and optimization of Systems.Antibes,FR(in French).D.R.C OX and D.V.H INKLEY(1986).Theoretical Statistics.Chapman and Hall,New York.D.R.C OX and H.D.M ILLER(1965).The Theory of Stochastic Processes.Wiley,New York.S.V.C ROWDER(1987).A simple method for studying run-length distributions of exponentially weighted moving average charts.Technometrics,vol.29,no4,pp.401-407.H.C S¨ORG¨O and L.H ORV´ATH(1988).Nonparametric methods for change point problems.In Handbook of Statistics(P.R.Krishnaiah,C.R.Rao,eds.),vol.7,Elsevier,New York,pp.403-425.R.B.D AVIES(1973).Asymptotic inference in stationary gaussian time series.Advances Applied Probability, vol.5,no3,pp.469-497.J.C.D ECKERT,M.N.D ESAI,J.J.D EYST and A.S.W ILLSKY(1977).F-8DFBW sensor failure identiﬁcation using analytical redundancy.IEEE Trans.Automatic Control,vol.AC-22,no5,pp.795-803.M.H.D E G ROOT(1970).Optimal Statistical Decisions.Series in Probability and Statistics,McGraw-Hill, New York.J.D ESHAYES and D.P ICARD(1979).Tests de ruptures dans un mod`e pte-Rendus de l’Acad´e mie des Sciences,vol.288,Ser.A,pp.563-566(in French).B IBLIOGRAPHY429 J.D ESHAYES and D.P ICARD(1983).Ruptures de Mod`e les en Statistique.Th`e ses d’Etat,Universit´e deParis-Sud,Orsay,France(in French).J.D ESHAYES and D.P ICARD(1986).Off-line statistical analysis of change-point models using non para-metric and likelihood methods.In Detection of Abrupt Changes in Signals and Dynamical Systems(M. Basseville,A.Benveniste,eds.).Lecture Notes in Control and Information Sciences,LNCIS77,Springer, New York,pp.103-168.B.D EVAUCHELLE-G ACH(1991).Diagnostic M´e canique des Fatigues sur les Structures Soumises`a des Vibrations en Ambiance de Travail.Th`e se de l’Universit´e Paris IX Dauphine(in French).B.D EVAUCHELLE-G ACH,M.B ASSEVILLE and A.B ENVENISTE(1991).Diagnosing mechanical changes in vibrating systems.Proc.SAFEPROCESS’91,Baden-Baden,FRG,pp.85-89.R.D I F RANCESCO(1990).Real-time speech segmentation using pitch and convexity jump models:applica-tion to variable rate speech coding.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-38,no5, pp.741-748.X.D ING and P.M.F RANK(1990).Fault detection via factorization approach.Systems and Control Letters, vol.14,pp.431-436.J.L.D OOB(1953).Stochastic Processes.Wiley,New York.V.D RAGALIN(1988).Asymptotic solutions in detecting a change in distribution under an unknown param-eter.Statistical Problems of Control,Issue83,Vilnius,pp.45-52.B.D UBUISSON(1990).Diagnostic et Reconnaissance des Formes.Trait´e des Nouvelles Technologies,S´e rie Diagnostic et Maintenance,Herm`e s,Paris(in French).A.J.D UNCAN(1986).Quality Control and Industrial Statistics,5th edition.Richard D.Irwin,Inc.,Home-wood,IL.J.D URBIN(1971).Boundary-crossing probabilities for the Brownian motion and Poisson processes and techniques for computing the power of the Kolmogorov-Smirnov test.Jal Applied Probability,vol.8,pp.431-453.J.D URBIN(1985).Theﬁrst passage density of the crossing of a continuous Gaussian process to a general boundary.Jal Applied Probability,vol.22,no1,pp.99-122.A.E MAMI-N AEINI,M.M.A KHTER and S.M.R OCK(1988).Effect of model uncertainty on failure detec-tion:the threshold selector.IEEE Trans.Automatic Control,vol.AC-33,no12,pp.1106-1115.J.D.E SARY,F.P ROSCHAN and D.W.W ALKUP(1967).Association of random variables with applications. Annals Mathematical Statistics,vol.38,pp.1466-1474.W.D.E WAN and K.W.K EMP(1960).Sampling inspection of continuous processes with no autocorrelation between successive results.Biometrika,vol.47,pp.263-280.G.F AVIER and A.S MOLDERS(1984).Adaptive smoother-predictors for tracking maneuvering targets.Proc. 23rd Conf.Decision and Control,Las Vegas,NV,pp.831-836.W.F ELLER(1966).An Introduction to Probability Theory and Its Applications,vol.2.Series in Probability and Mathematical Statistics,Wiley,New York.R.A.F ISHER(1925).Theory of statistical estimation.Proc.Cambridge Philosophical Society,vol.22, pp.700-725.M.F ISHMAN(1988).Optimization of the algorithm for the detection of a disorder,based on the statistic of exponential smoothing.In Statistical Problems of Control,Issue83,Vilnius,pp.146-151.R.F LETCHER(1980).Practical Methods of Optimization,2volumes.Wiley,New York.P.M.F RANK(1990).Fault diagnosis in dynamic systems using analytical and knowledge based redundancy -A survey and new results.Automatica,vol.26,pp.459-474.430B IBLIOGRAPHY P.M.F RANK(1991).Enhancement of robustness in observer-based fault detection.Proc.SAFEPRO-CESS’91,Baden-Baden,FRG,pp.275-287.P.M.F RANK and J.W¨UNNENBERG(1989).Robust fault diagnosis using unknown input observer schemes. In Fault Diagnosis in Dynamic Systems-Theory and Application(R.Patton,P.Frank,R.Clark,eds.). International Series in Systems and Control Engineering,Prentice Hall International,London,UK,pp.47-98.K.F UKUNAGA(1990).Introduction to Statistical Pattern Recognition,2d ed.Academic Press,New York. S.I.G ASS(1958).Linear Programming:Methods and Applications.McGraw Hill,New York.W.G E and C.Z.F ANG(1989).Extended robust observation approach for failure isolation.Int.Jal Control, vol.49,no5,pp.1537-1553.W.G ERSCH(1986).Two applications of parametric time series modeling methods.In Mechanical Signature Analysis-Theory and Applications(S.Braun,ed.),chap.10.Academic Press,London.J.J.G ERTLER(1988).Survey of model-based failure detection and isolation in complex plants.IEEE Control Systems Magazine,vol.8,no6,pp.3-11.J.J.G ERTLER(1991).Analytical redundancy methods in fault detection and isolation.Proc.SAFEPRO-CESS’91,Baden-Baden,FRG,pp.9-22.B.K.G HOSH(1970).Sequential Tests of Statistical Hypotheses.Addison-Wesley,Cambridge,MA.I.N.G IBRA(1975).Recent developments in control charts techniques.Jal Quality Technology,vol.7, pp.183-192.J.P.G ILMORE and R.A.M C K ERN(1972).A redundant strapdown inertial reference unit(SIRU).Jal Space-craft,vol.9,pp.39-47.M.A.G IRSHICK and H.R UBIN(1952).A Bayes approach to a quality control model.Annals Mathematical Statistics,vol.23,pp.114-125.A.L.G OEL and S.M.W U(1971).Determination of the ARL and a contour nomogram for CUSUM charts to control normal mean.Technometrics,vol.13,no2,pp.221-230.P.L.G OLDSMITH and H.W HITFIELD(1961).Average run lengths in cumulative chart quality control schemes.Technometrics,vol.3,pp.11-20.G.C.G OODWIN and K.S.S IN(1984).Adaptive Filtering,Prediction and rmation and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.R.M.G RAY and L.D.D AVISSON(1986).Random Processes:a Mathematical Approach for Engineers. Information and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.C.G UEGUEN and L.L.S CHARF(1980).Exact maximum likelihood identiﬁcation for ARMA models:a signal processing perspective.Proc.1st EUSIPCO,Lausanne.D.E.G USTAFSON, A.S.W ILLSKY,J.Y.W ANG,M.C.L ANCASTER and J.H.T RIEBWASSER(1978). ECG/VCG rhythm diagnosis using statistical signal analysis.Part I:Identiﬁcation of persistent rhythms. Part II:Identiﬁcation of transient rhythms.IEEE Trans.Biomedical Engineering,vol.BME-25,pp.344-353 and353-361.F.G USTAFSSON(1991).Optimal segmentation of linear regression parameters.Proc.IFAC/IFORS Symp. Identiﬁcation and System Parameter Estimation,Budapest,pp.225-229.T.H¨AGGLUND(1983).New Estimation Techniques for Adaptive Control.Ph.D.Thesis,Lund Institute of Technology,Lund,Sweden.T.H¨AGGLUND(1984).Adaptive control of systems subject to large parameter changes.Proc.IFAC9th World Congress,Budapest.B IBLIOGRAPHY431 P.H ALL and C.C.H EYDE(1980).Martingale Limit Theory and its Application.Probability and Mathemat-ical Statistics,a Series of Monographs and Textbooks,Academic Press,New York.W.J.H ALL,R.A.W IJSMAN and J.K.G HOSH(1965).The relationship between sufﬁciency and invariance with applications in sequential analysis.Ann.Math.Statist.,vol.36,pp.576-614.E.J.H ANNAN and M.D EISTLER(1988).The Statistical Theory of Linear Systems.Series in Probability and Mathematical Statistics,Wiley,New York.J.D.H EALY(1987).A note on multivariate CuSum procedures.Technometrics,vol.29,pp.402-412.D.M.H IMMELBLAU(1970).Process Analysis by Statistical Methods.Wiley,New York.D.M.H IMMELBLAU(1978).Fault Detection and Diagnosis in Chemical and Petrochemical Processes. Chemical Engineering Monographs,vol.8,Elsevier,Amsterdam.W.G.S.H INES(1976a).A simple monitor of a system with sudden parameter changes.IEEE r-mation Theory,vol.IT-22,no2,pp.210-216.W.G.S.H INES(1976b).Improving a simple monitor of a system with sudden parameter changes.IEEE rmation Theory,vol.IT-22,no4,pp.496-499.D.V.H INKLEY(1969).Inference about the intersection in two-phase regression.Biometrika,vol.56,no3, pp.495-504.D.V.H INKLEY(1970).Inference about the change point in a sequence of random variables.Biometrika, vol.57,no1,pp.1-17.D.V.H INKLEY(1971).Inference about the change point from cumulative sum-tests.Biometrika,vol.58, no3,pp.509-523.D.V.H INKLEY(1971).Inference in two-phase regression.Jal American Statistical Association,vol.66, no336,pp.736-743.J.R.H UDDLE(1983).Inertial navigation system error-model considerations in Kalmanﬁltering applica-tions.In Control and Dynamic Systems(C.T.Leondes,ed.),Academic Press,New York,pp.293-339.J.S.H UNTER(1986).The exponentially weighted moving average.Jal Quality Technology,vol.18,pp.203-210.I.A.I BRAGIMOV and R.Z.K HASMINSKII(1981).Statistical Estimation-Asymptotic Theory.Applications of Mathematics Series,vol.16.Springer,New York.R.I SERMANN(1984).Process fault detection based on modeling and estimation methods-A survey.Auto-matica,vol.20,pp.387-404.N.I SHII,A.I WATA and N.S UZUMURA(1979).Segmentation of nonstationary time series.Int.Jal Systems Sciences,vol.10,pp.883-894.J.E.J ACKSON and R.A.B RADLEY(1961).Sequential and tests.Annals Mathematical Statistics, vol.32,pp.1063-1077.B.J AMES,K.L.J AMES and D.S IEGMUND(1988).Conditional boundary crossing probabilities with appli-cations to change-point problems.Annals Probability,vol.16,pp.825-839.M.K.J EERAGE(1990).Reliability analysis of fault-tolerant IMU architectures with redundant inertial sen-sors.IEEE Trans.Aerospace and Electronic Systems,vol.AES-5,no.7,pp.23-27.N.L.J OHNSON(1961).A simple theoretical approach to cumulative sum control charts.Jal American Sta-tistical Association,vol.56,pp.835-840.N.L.J OHNSON and F.C.L EONE(1962).Cumulative sum control charts:mathematical principles applied to their construction and use.Parts I,II,III.Industrial Quality Control,vol.18,pp.15-21;vol.19,pp.29-36; vol.20,pp.22-28.432B IBLIOGRAPHY R.A.J OHNSON and M.B AGSHAW(1974).The effect of serial correlation on the performance of CUSUM tests-Part I.Technometrics,vol.16,no.1,pp.103-112.H.L.J ONES(1973).Failure Detection in Linear Systems.Ph.D.Thesis,Dept.Aeronautics and Astronautics, MIT,Cambridge,MA.R.H.J ONES,D.H.C ROWELL and L.E.K APUNIAI(1970).Change detection model for serially correlated multivariate data.Biometrics,vol.26,no2,pp.269-280.M.J URGUTIS(1984).Comparison of the statistical properties of the estimates of the change times in an autoregressive process.In Statistical Problems of Control,Issue65,Vilnius,pp.234-243(in Russian).T.K AILATH(1980).Linear rmation and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.L.V.K ANTOROVICH and V.I.K RILOV(1958).Approximate Methods of Higher Analysis.Interscience,New York.S.K ARLIN and H.M.T AYLOR(1975).A First Course in Stochastic Processes,2d ed.Academic Press,New York.S.K ARLIN and H.M.T AYLOR(1981).A Second Course in Stochastic Processes.Academic Press,New York.D.K AZAKOS and P.P APANTONI-K AZAKOS(1980).Spectral distance measures between gaussian pro-cesses.IEEE Trans.Automatic Control,vol.AC-25,no5,pp.950-959.K.W.K EMP(1958).Formula for calculating the operating characteristic and average sample number of some sequential tests.Jal Royal Statistical Society,vol.B-20,no2,pp.379-386.K.W.K EMP(1961).The average run length of the cumulative sum chart when a V-mask is used.Jal Royal Statistical Society,vol.B-23,pp.149-153.K.W.K EMP(1967a).Formal expressions which can be used for the determination of operating character-istics and average sample number of a simple sequential test.Jal Royal Statistical Society,vol.B-29,no2, pp.248-262.K.W.K EMP(1967b).A simple procedure for determining upper and lower limits for the average sample run length of a cumulative sum scheme.Jal Royal Statistical Society,vol.B-29,no2,pp.263-265.D.P.K ENNEDY(1976).Some martingales related to cumulative sum tests and single server queues.Stochas-tic Processes and Appl.,vol.4,pp.261-269.T.H.K ERR(1980).Statistical analysis of two-ellipsoid overlap test for real time failure detection.IEEE Trans.Automatic Control,vol.AC-25,no4,pp.762-772.T.H.K ERR(1982).False alarm and correct detection probabilities over a time interval for restricted classes of failure detection algorithms.IEEE rmation Theory,vol.IT-24,pp.619-631.T.H.K ERR(1987).Decentralizedﬁltering and redundancy management for multisensor navigation.IEEE Trans.Aerospace and Electronic systems,vol.AES-23,pp.83-119.Minor corrections on p.412and p.599 (May and July issues,respectively).R.A.K HAN(1978).Wald’s approximations to the average run length in cusum procedures.Jal Statistical Planning and Inference,vol.2,no1,pp.63-77.R.A.K HAN(1979).Someﬁrst passage problems related to cusum procedures.Stochastic Processes and Applications,vol.9,no2,pp.207-215.R.A.K HAN(1981).A note on Page’s two-sided cumulative sum procedures.Biometrika,vol.68,no3, pp.717-719.B IBLIOGRAPHY433 V.K IREICHIKOV,V.M ANGUSHEV and I.N IKIFOROV(1990).Investigation and application of CUSUM algorithms to monitoring of sensors.In Statistical Problems of Control,Issue89,Vilnius,pp.124-130(in Russian).G.K ITAGAWA and W.G ERSCH(1985).A smoothness prior time-varying AR coefﬁcient modeling of non-stationary covariance time series.IEEE Trans.Automatic Control,vol.AC-30,no1,pp.48-56.N.K LIGIENE(1980).Probabilities of deviations of the change point estimate in statistical models.In Sta-tistical Problems of Control,Issue83,Vilnius,pp.80-86(in Russian).N.K LIGIENE and L.T ELKSNYS(1983).Methods of detecting instants of change of random process prop-erties.Automation and Remote Control,vol.44,no10,Part II,pp.1241-1283.J.K ORN,S.W.G ULLY and A.S.W ILLSKY(1982).Application of the generalized likelihood ratio algorithm to maneuver detection and estimation.Proc.American Control Conf.,Arlington,V A,pp.792-798.P.R.K RISHNAIAH and B.Q.M IAO(1988).Review about estimation of change points.In Handbook of Statistics(P.R.Krishnaiah,C.R.Rao,eds.),vol.7,Elsevier,New York,pp.375-402.P.K UDVA,N.V ISWANADHAM and A.R AMAKRISHNAN(1980).Observers for linear systems with unknown inputs.IEEE Trans.Automatic Control,vol.AC-25,no1,pp.113-115.S.K ULLBACK(1959).Information Theory and Statistics.Wiley,New York(also Dover,New York,1968). K.K UMAMARU,S.S AGARA and T.S¨ODERSTR¨OM(1989).Some statistical methods for fault diagnosis for dynamical systems.In Fault Diagnosis in Dynamic Systems-Theory and Application(R.Patton,P.Frank,R. Clark,eds.).International Series in Systems and Control Engineering,Prentice Hall International,London, UK,pp.439-476.A.K USHNIR,I.N IKIFOROV and I.S AVIN(1983).Statistical adaptive algorithms for automatic detection of seismic signals-Part I:One-dimensional case.In Earthquake Prediction and the Study of the Earth Structure,Naouka,Moscow(Computational Seismology,vol.15),pp.154-159(in Russian).L.L ADELLI(1990).Diffusion approximation for a pseudo-likelihood test process with application to de-tection of change in stochastic system.Stochastics and Stochastics Reports,vol.32,pp.1-25.T.L.L A¨I(1974).Control charts based on weighted sums.Annals Statistics,vol.2,no1,pp.134-147.T.L.L A¨I(1981).Asymptotic optimality of invariant sequential probability ratio tests.Annals Statistics, vol.9,no2,pp.318-333.D.G.L AINIOTIS(1971).Joint detection,estimation,and system identiﬁrmation and Control, vol.19,pp.75-92.M.R.L EADBETTER,G.L INDGREN and H.R OOTZEN(1983).Extremes and Related Properties of Random Sequences and Processes.Series in Statistics,Springer,New York.L.L E C AM(1960).Locally asymptotically normal families of distributions.Univ.California Publications in Statistics,vol.3,pp.37-98.L.L E C AM(1986).Asymptotic Methods in Statistical Decision Theory.Series in Statistics,Springer,New York.E.L.L EHMANN(1986).Testing Statistical Hypotheses,2d ed.Wiley,New York.J.P.L EHOCZKY(1977).Formulas for stopped diffusion processes with stopping times based on the maxi-mum.Annals Probability,vol.5,no4,pp.601-607.H.R.L ERCHE(1980).Boundary Crossing of Brownian Motion.Lecture Notes in Statistics,vol.40,Springer, New York.L.L JUNG(1987).System Identiﬁcation-Theory for the rmation and System Sciences Series, Prentice Hall,Englewood Cliffs,NJ.。

gmm

1.2 Single Equation Linear GMM
Consider the linear regression model yt = z0 t δ 0 + εt , t = 1, . . . , n (1.1)
where zt is an L × 1 vector of explanatory variables, δ 0 is a vector of unknown coeﬃcients and εt is a random error term. The model (1.1) allows for the possibility that some or all of the elements of zt may be correlated with the error term εt , i.e., E [ztk εt ] 6= 0 for some k. If E [ztk εi ] 6= 0 then ztk is called an endogenous variable . It is well known that if zt contains endogenous variables then the least squares estimator of δ 0 in (1.1) is biased and inconsistent. Associated with the model (1.1), it is assumed that there exists a K × 1 vector of instrumental variables xt which may contain some or all of the elements of zt . Let wt represent the vector of unique and non-constant elements of {yt , zt , xt }. It is assumed that {wt } is a stationary and ergodic stochastic process. The instrumental variables xt satisfy the set of K orthogonality conditions E [gt (wt , δ 0 )] = E [xt εt ] = E [xt (yt − z0 (1.2) t δ 0 )] = 0 where gt (wt , δ 0 ) = xt εt = xt (yt − z0 t δ 0 ). Expanding (1.2), gives the relation Σxy = Σxz δ 0 where Σxy = E [xt yt ] and Σxz = E [xt z0 t ]. For identiﬁcation of δ 0 , it is required that the K × L matrix E [xt z0 t ] = Σxz be of full rank L. This rank condition ensures that δ 0 is the unique solution to (1.2). Note, if K = L, then Σxz is invertible and δ 0 may be determined using

Volatility in Financial Markets Stochastic Models and Empirical Results

a r X i v :c o n d -m a t /0202527v 1 [c o n d -m a t .s t a t -m e c h ] 28 F eb 2002Volatility in Financial Markets:Stochastic Models and Empirical Results Salvatore Miccich`e ,Giovanni Bonanno,Fabrizio Lillo and Rosario N.Mantegna 1Istituto Nazionale per la Fisica della Materia,Unit`a di Palermo and Dipartimento di Fisica e Tecnologie Relative,Universit`a di Palermo,Viale delle Scienze,I-90128,Palermo,Italia 1Introduction Volatility of ﬁnancial time series is a key variable in the modeling of ﬁnancial markets.It controls all the risk measures associated with the dynamics of price of a ﬁnancial asset.It also aﬀects the rational price of derivative products.In this paper we consider some stochastic volatility models proposed in the ﬁnan-cial literature by investigating their ability in modeling statistical properties detected in empirical data.Speciﬁcally,we investigate the probability density function (pdf)of historical volatility for 100highly capitalized stocks traded in the US equity markets.We observe that widespread volatility models such as the Hull and White model [1]and the lognormal model fail in describingthe volatility pdf when we ask the model to describe both low and high values of volatility.Our results show that a lognormal pdf better describes low values of volatility whereas the Hull and White pdf gives a better approximation of the empirical pdf for large values.2Volatility modelsThe volatilityσof aﬁnancial asset is a statistical quantity which needs to be determined starting from market information[2].It is the standard devi-ation of asset return(or,almost equivalently,of logarithm price changes of the asset).Diﬀerent methodologies are used to infer volatility estimation from market data ranging from a direct calculation from past return data(historical volatility)to the computation of the volatility implied in the determination of an option price computed using the Black and Scholes formula[3]or some vari-ant of it.There is a large empirical evidence that volatility is itself a stochastic process.In the present study we aim at comparing the theoretical predictions for the pdf of the volatilityσobtained with two diﬀerent stochastic volatil-ity models with empirical observations obtained for the100most capitalized stocks traded in US equity markets(mostly the New York Stock Exchange and the NASDAQ).Theﬁrst model we consider is the Hull and White model [1].In this model,the variance rate v≡σ2is described by the Ito’s equationdv=a(b−v)dt+ξvαdz v(1) where a and b are parameters controlling the mean reverting nature of the stochastic process,ξis controlling its diﬀusive aspects and z v is a Wiener pro-cess.The stochastic process is reverting at a level b at a rate a.The exponent αhas been set to1in the investigation of Hull and White[1]and to1/2in the investigation of Heston[4].In the present study we investigate the Hull and White model withα=1.The Ito’s equation of this model for the volatilityσisdσ=12ξ2σ2 dt+ξσ2dzσ (2)This equation has been obtained starting from Eq.(1)and using Ito’s lemma. The Hull and White model has associated a stationary pdf of the volatility which has the formP(σ)=2(ba/ξ2)1+a/ξ2σ2a/ξ2+3(3)0246810σn 00.511.5P (σn )024681000.511.5P (σn )0246810σn 10−510−410−310−210−1100101024681010−510−410−310−210−1100101(a)(c)Fig.1.Best ﬁts of the empirical pdf of normalized volatility obtained by investigating 100stocks traded in US equity markets during the time period January 1995-December 1998.In panel (a)and (c)the y axis is linear whereas in panel in (b)and (d)the data are shown in a semilogarithmic plot.In all panels,the solid lines are the best ﬁts whereas the histogram and solid circles are empirical data.In panels (a)and (b)we show the best ﬁttings obtained with a lognormal pdf of mean value 0.97and variance equals to 0.19and in panels (c)and (d)we show the best ﬁttings obtained with the Hull and White pdf of Eq.(3).In this last case the ﬁttings parameters are 2a/ξ2+3=3.79and ba/ξ2=0.91.This pdf has a power-law tail for large values of σ.A power-law tail in the empirical volatility pdf has been observed in Ref.[5]for large values of the volatility.Another model is the lognormal model of volatility [6,7].An Ito’s stochastic diﬀerential equation associated with a lognormal pdf isdσ=a (b −ln σ)dt +ξσ1/2dz σ(4)where a ,b and ξare control parameters of the model.The two models are characterized by quite diﬀerent pdfs especially for large values of the volatility where the Hull and White pdf shows a power-law behavior.The present study aims at detecting the regions of validity of these two models in empirical data.n 00.51P (σn )Fig.2.Empirical pdf of normalized volatility compared with the pdf predicted by the stochastic model of Eq.(5).In the ﬁgure the y axis is linear whereas in the inset the data are shown in a log-log plot.Solid lines represent the best ﬁt whereas the histogram and solid circles are empirical data.The ﬁttings parameters of Eq.(6)are a/ξ2+2=6.27and ba/ξ2=4.45.3Empirical ResultsWith this goal in mind,we investigate the statistical properties of volatility for the 100most capitalized stocks traded in US equity markets during a 4year time period.The empirical data are taken from the trade and quote (TAQ)database,maintained by the New York Stock Exchange (NYSE).In particular,our data cover the whole period ranging from January 1995to December 1998(1011trading days).This database contains all transactions occurred for each stock traded in the US equity markets.The capitalization considered is the one recorded on August 31,1998.For each stock and for each trading day we consider the time series of stock price recorded transaction by transaction.Since transactions for diﬀerent stocks do not happen simultaneously,we divide each trading day (lasting 6h 30′)into 12intervals of 1950seconds each.In correspondence to each interval,we deﬁne 12(intraday)stock’s prices proxies S (k )–with k =1,···,12deﬁned as the transaction price detected nearest to end of the interval (this a one possible way to deal with high-frequency ﬁnancial data [8]).We choose 12intraday intervals since this value ensures that at least one transaction is in average observed in each interval for all consideredstocks in the present study.For each stock we can thus compute a historical√daily volatility asσ(t)=exp(−ba/ξ2σ)Γ(1+a/ξ2)can be written in terms of a/ξ2and ba/ξ2.In Fig.2we show the bestﬁt with the pdf of Eq.(6).The agreement with empirical data is rather good in a volatility range fromσn=0.5toσn=10.In summary,we report on a comparison of two widespread theoretical models of volatility with empirical data obtained by collecting together the volatility of100most capitalized stocks traded in US equity markets.The comparison is focused on the shape of the asymptotic pdf of volatility.Two widespread models(lognormal and Hull and White)fail in describing the pdf over a rel-atively wide volatility range.We show that the model of Eq.(5)improves the overall description of the pdf especially for values of normalized volatility σn>0.5.Further research attempts are needed to select the most appropri-ate Ito’s model able to describe volatility both under the aspects of the pdf and under the dynamics aspects of the nature and form of volatility auto-correlation function.Indeed,there is a growing evidence that the volatility autocorrelation function is long-range correlated[8]and this key aspect is not taken into account in most of the widespread models of volatility(as the ones considered in the present study)which are typically characterized only by short-range time memory.AcknowledgementsThe authors thank INFM,MIUR and ASI forﬁnancial support.This work is part of the FRA-INFM project’Volatility inﬁnancial markets’. References[1]J.Hull,A.White,Journal of Finance XLII(1987)281.[2]J.C.Hull Options,Futures and Other Derivatives Prentice Hall,Inc.(1997)[3] F.Black,M.Scholes,Journal of Political Economy81(1973)637.[4]S.L.Heston,Review of Financial Studies6(1993)327.[5]Y.Liu,P.Gopikrishnan,P.Cizeau,M.Meyer,C.-K.Peng,H.E.Stanley,Phys.Rev.D60(1999)1390.[6]P.Cizeau,Y.Liu,M.Meyer,C.-K.Peng,and H.E.Stanley,Physica A245(1997)441.[7]M.Pasquini,M.Serva,Physica A269(1999)140.[8]M.M.Dacorogna,R.Gencay,U.A.M¨u ller,R.B.Olsen,O.V.Pictet,AnIntroduction to High-Frequency Finance,Academic Press(2001)。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。