如题,如果我想重覆一段R code若干次(如100000次),要如何把我写的code放进同一个loop里面呢,让结果能够自动化的存进一个表内?
下面是code的内容:
## Fixed quantities
# One-step transition probability matrix P
#_______________________________________________________________________________
P=matrix(c(0.8,0.1,0.1,
0.3,0.65,0.05,
0.3,0.05,0.65), nrow=3, ncol=3, byrow=T)
P
# Boundaries
#_______________________________________________________________________________
B1_state1=qgamma(0.3, shape = 1, rate = 0.1, lower.tail = TRUE)
B2_state1=qgamma(0.7, shape = 1, rate = 0.1, lower.tail = TRUE)
B1_state2=qgamma(0.3, shape = 0.5, rate = 0.1, lower.tail = TRUE)
B2_state2=qgamma(0.7, shape = 0.5, rate = 0.1, lower.tail = TRUE)
B1_state3=qgamma(0.3, shape = 0.9, rate = 0.06, lower.tail = TRUE)
B2_state3=qgamma(0.7, shape = 0.9, rate = 0.06, lower.tail = TRUE)
# Premium loading c
#_______________________________________________________________________________
C=c(120,140,160,180,200)/100
## Vary quantities
# Aggregate claims distribution
#_______________________________________________________________________________
aggregate_claim_in_state_1=rgamma(50, shape = 1, rate = 0.1) #state 1
aggregate_claim_in_state_2=rgamma(50, shape = 0.5, rate = 0.1) #state 2
aggregate_claim_in_state_3=rgamma(50, shape = 0.9, rate = 0.06) #state3
# Economic states
#_______________________________________________________________________________
simMarkov = function(P, lens=50){
n=nrow(P)
result = numeric(lens)
result[1]=1
for (i in 2:lens){
result[i]=sample(1:n, 1, prob=P[result[i-1],])
}
result
}
Econ_state = simMarkov(P)
# Lower bound & upper bound
#_______________________________________________________________________________
lower.bound=vector(length = 50)
for (i in 1:50) {
if(Econ_state[i]==1){
lower.bound[i]=B1_state1}
else if(Econ_state[i]==2) {
lower.bound[i]=B1_state2}
else if(Econ_state[i]==3) {
lower.bound[i]=B1_state3}
}
upper.bound=vector(length = 50)
for (i in 1:50) {
if(Econ_state[i]==1){
upper.bound[i]=B2_state1}
else if(Econ_state[i]==2) {
upper.bound[i]=B2_state2}
else if(Econ_state[i]==3) {
upper.bound[i]=B2_state3}
}
# Aggregate claims
#_______________________________________________________________________________
s=vector(length = 50)
for (k in 1:50) {
if(Econ_state[k] == 1){
s[k]=aggregate_claim_in_state_1[k]
} else if(Econ_state[k] == 2){
s[k]=aggregate_claim_in_state_2[k]
} else if(Econ_state[k] == 3){
s[k]=aggregate_claim_in_state_3[k]
}
}
# Expected Aggregated claims
#_______________________________________________________________________________
mean_s=vector(length = 50)
for (i in 1:50) {
expected_aggregate_claim=
c(aggregate_claim_in_state_1[i],aggregate_claim_in_state_2[i],
aggregate_claim_in_state_3[i])
mean_s[i]=mean(expected_aggregate_claim)
}
# Premium loading
#_______________________________________________________________________________
c1=vector(length = 51)
for (i in 2:51) {
c1[1]=C[1]
if(s[i-1]<lower.bound[i-1]){
c1[i]=c1[i-1]-0.2
}else if(lower.bound[i-1]<=s[i-1] && s[i-1]<upper.bound[i-1]){
c1[i]=c1[i-1]
}else if(s[i-1]>=upper.bound[i-1] && c1[i-1]>=2){
c1[i]=C[5]
}else if (s[i-1]<lower.bound[i-1] && c1[i-1]<=1.2) {
c1[i]=C[1]
}else if (s[i-1]>=upper.bound[i-1]) {
c1[i]=c1[i-1]+0.2
}
}
# Surplus level
#_______________________________________________________________________________
U=matrix(,nrow = 50, ncol = 201)
initial.u=0:200
for (n in 2:50) {
for (u in 1:201) {
u0=initial.u[u]
U[1,u]=u0+c1[2]*mean_s[1]-s[1]
U[n,u]=U[n-1,u]+c1[n+1]*mean_s[n]-s[n]
}
}
# Ruin Probability
#_______________________________________________________________________________
Ruin1=vector(length = 201)
for (j in 1:201) {
Ruin1[j]=sum(U[,j]<0)/50
}
Ruin1