Codes/Figures                        Données                         Compléments                     Erratum





Comp7. Simulation



# 7.4.2 Simulation d'une intervention

impulse.plot=function(x,y,couleur= NA, lwd =2){
# dessine les impulsions de valeur y en x
# x abscisse, vecteur
# y ordonnée, vecteur de mm dimension
# si y[i] non nul trait vertical de (x[i],0) à (x[i], y[i])
np = length(x)
num = 1:np
i.n0 =  num[y != 0]
if (is.na(couleur) == TRUE)  { couleur = "red"}
for (j in i.n0)
{
# segment vertical (1) abscisse, ordonnée du  départ, (2) abscisse ordonnée de l'arrivée
lines(c(x[j],x[j]), c(0, y[j]),col=couleur, lwd=2)
}
}

imp.fun= function(temps,t0)
{a=rep(0,length(temps)); a[t0]=1 ;a }
t0 = 3
dates=1:20; ldat=length(dates)
y.imp = imp.fun(dates,t0)
delta = .8; omega = 4 ; vecauto= c(1,-delta)
AR =array( vecauto,c(length(vecauto),1,1))
vecma = 1; MA=array(vecma,c(length(vecma),1,1))

# rappel (chapitre 7)
# pour obtenir y1
require("dse")
mod1 = ARMA(A=AR,B=MA)
y1 = as.vector(simulate(mod1,y0=0,noise=as.matrix(omega*y.imp),
        sampleT=ldat)$output)

# rappel (chapitre 7)
# pour obtenir y2
ech.fun= function(temps,t0)
{a=rep(0,length(temps)); a[temps >= t0]=1 ;a }
y.ech = ech.fun(dates,t0)
y2 = as.vector(simulate(mod1,y0=0,noise=as.matrix(omega*y.ech),
          sampleT=ldat)$output)



mar.int1 = c(1,1,0,.5)  # graphe du bas
mar.int2 = c(0,1,.5,.5)  # graphe du haut
# h.in = hauteur en pouce
# mb.pouce = marge du bas en pouce
# hauteur totale en pouce
h.in=4
# moins les marges
# en dehors marge du bas et du bas , il reste en tout
h.reste.t = h.in - mar.int1[1] - mar.int2[3]
# hauteurs de bas en haut
hg1 = 0.6*h.reste.t
hg2 = 0.4*h.reste.t
# haut du graphe du bas
pr1 = (hg1+mar.int1[1])/h.in
# mgp The margin line (in mex units) for the axis title, axis labels and axis line.
op = par(oma=c(0,0,0,0),mgp=c(2.7,.7,0),cex=.9,cex.lab=1,cex.axis=.9 )
op1 = par(fig= c(0,1,pr1,1), mai=mar.int2)  # mai remplace mar
plot(dates,y1,type="n", xaxt="n",ylab="réponse",xlab="")
impulse.plot(dates,y1)
text(11,3,labels="Amortissement progressif d'une impulsion")
# chronogramme bas
par(new=T)
op2 = par(fig= c(0,1,0,pr1),mai=mar.int1)                                                                                       
plot(dates,y2,type="n",xaxt="n", ylab="réponse",xlab="temps")
impulse.plot(dates,y2)
axis(1,at= c(1,3,seq(5,20,by=5)), lab=c("1","t0","5","10","15","20"))
text(1,16,labels="Réponse progressive à un échelon",pos=4)
par(op2)
par(op1)
par(op)

Fig7-1