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
84
85
86
87
88
89
90
91
92
93
94
|
interpweights <- function(w, v1, v2){
#Given L=(w,v1), compute neww such that newL=(new,v2)=L in distribution
cumw <- cumsum(w)
neww <- splinefun(v1,cumw,method= "monoH.FC")(v2,deriv=1)
#neww <- diff(newcumw)
interpweights <- neww/sum(neww)
return(interpweights)
}
adjust_scenario <- function(scenario, epsilon){
1-(1-scenario)^(1/(1+epsilon))
}
adjust_weights <- function(weights, scenario, epsilon){
interpweights(weights,scenario,adjust_scenario(scenario,epsilon))
}
obj <- function(epsilon, vecpv, prob, support, cte){
newprob <- adjust_weights(prob, support, epsilon)
return( 1 - crossprod(newprob, vecpv) - cte)
}
optimize <- function(min, max, vecpv, prob, support, cte){
mid <- (min + max)/2
objective <- obj(mid, vecpv, prob, support, cte)
while( abs(objective)>1e-6){
if(objective>0){
min <- mid
}else{
max <- mid
}
mid <- (min+max)/2
objective <- obj(mid, vecpv, prob, support, cte)
}
return( mid )
}
interpweightsadjust <- function(w, v1, v2, vecpv){
interpweightsadjust <- interpweights(w, v1, v2)
epsilon <- optimize(-0.5, 0.5, vecpv, interpweightsadjust, v2, 1)
return( adjust_weights(interpweightsadjust, v2, epsilon) )
}
transformweightslike <- function(p1, v1, p2, v2, p, v){
cump2 <- cumsum(p2)
cump1 <- cumsum(p1)
P1 <- splinefun(v1,cump1,method= "monoH.FC")
dP1 <- function(x){P1(x,deriv=1)}
pomme <- interpweights(p2,v2,v)
pomme <- cumsum(pomme)
r <- rep(0,length(pomme))
for(i in 1:length(pomme)){
r[i] <- inverse(P1,dP1,pomme[i])
}
return(r)
}
clipw <- function(x){
write(x,file="clipboard",sep="\n")
}
clipr <- function(){
scan(file="clipboard")
}
sclipr <- function(){
scan(file="clipboard",what="character")
}
inverse <- function(f,Df,x, x0=x){
#inverse a function by the newton's method.
x1 <- x0-f(x0)/(Df(x0)-1)
counter <- 0
while(abs(x1-x0)>1e-6&&counter<500){
x0 <- x1
x1 <- x0-(f(x0)-x)/Df(x0)
counter <- counter+1
}
return(x1)
}
setwd("W:/corpCDOS/NextGen")
lcdx10 <- read.table("lcdx10",sep=",",header=T)
lcdx12 <- read.table("lcdx12",sep=",",header=T)
hy10_5y <- read.table("hy10_5y",sep=",",header=T)
hy10_7y <- read.table("hy10_7y",sep=",",header=T)
lcdx10_pv <- read.table("lcdx10_pv",sep=",",header=T)
lcdx12_pv <- read.table("lcdx12_pv",sep=",",header=T)
hy10_5y_pv <- read.table("hy10_5y_pv",sep=",",header=T)
hy10_7y_pv <- read.table("hy10_7y_pv",sep=",",header=T)
clipw(interpweightsadjust(lcdx10$dist,lcdx10$support,lcdx10_pv$support,lcdx10_pv$pv))
clipw(interpweightsadjust(hy10_5y$dist,hy10_5y$support,hy10_5y_pv$support,hy10_5y_pv$pv))
clipw(interpweightsadjust(hy10_7y$dist,hy10_7y$support,hy10_7y_pv$support,hy10_7y_pv$pv))
clipw(interpweightsadjust(lcdx12$dist,lcdx12$support,lcdx12_pv$support,lcdx12_pv$pv))
|
