forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
19.R
83 lines (68 loc) · 2.09 KB
/
19.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
logistic.map <- function(x,r) {
return(4*r*x*(1-x))
}
logistic.map.ts <- function(timelength,r,initial.cond=NULL) {
x <-vector(mode="numeric",length=timelength)
if(is.null(initial.cond)) {
x[1] <-runif(1)
} else {
x[1] <-initial.cond
}
for (t in 2:timelength) {
x[t] = logistic.map(x[t-1],r)
}
return(x)
}
logistic.noisy.ts <- function(timelength,r,initial.cond=NULL,noise.sd=0.1) {
x <- logistic.map.ts(timelength,r,initial.cond)
return(x+rnorm(timelength,0,noise.sd))
}
mean.logistic <- function(r,n=1e3,s=10) {
mean(replicate(s,mean(logistic.map.ts(n,r))))
}
mean.logistic.plottable <- function(r,n=1e3,s=10) {
sapply(r,mean.logistic,n=n,s=s)
}
# curve(mean.logistic.plottable)
var.logistic <- function(r,n=1e3,s=10) {
mean(replicate(s,var(logistic.map.ts(n,r))))
}
var.logistic.plottable <- function(r,n=1e3,s=10) {
sapply(r,var.logistic,n=n,s=s)
}
# curve(var.logistic.plottable,from=0,to=1,xlab="r",ylab="Variance of X")
msm.logistic.est <- function(x,S=10) {
T <- length(x)
mom <- vector(length=2)
mom[1] <- mean(x)
mom[2] <- var(x)
moment.discrep <- function(r) {
sims <- replicate(S,logistic.map.ts(T,r))
momp <- vector(length=2)
momp[1] <- mean(apply(sims,2,mean))
momp[2] <- mean(apply(sims,2,var))
return(sum((mom-momp)^2))
}
return(optimize(f=moment.discrep,lower=0,upper=1)$minimum)
}
ar.from.logmap <- function(r,T=1000,reps=10,order=2,coef=1) {
mean(replicate(reps,ar(logistic.noisy.ts(T,r),aic=FALSE,order.max=order)$ar[coef]))
}
ar.logmap.plottable <- function(r,...) {sapply(r,ar.from.logmap,...)}
logistic.map.II <- function(y,order=2,S=10) {
T <- length(y)
ar.fit <- function(x) {
return(ar(x,aic=FALSE,order.max=order)$ar)
}
beta.data <- ar.fit(y)
beta.discrep <- function(r) {
beta.S <- mean(replicate(S,ar.fit(logistic.noisy.ts(T,r))))
return(sum((beta.data - beta.S)^2))
}
return(optimize(beta.discrep,lower=0.75,upper=1))
}
mse.logistic.II <- function(T,r=0.9,reps=30,order=2,S=10) {
II <- replicate(reps,logistic.map.II(logistic.noisy.ts(T,r),order=order,S=S)$minimum)
II.error = II - r # Uses duplication
return(mean(II.error^2))
}