# Forecast 2 year ahead rstan

I have done my work on forecast 1 year ahead for run-off triangle data but could not work out forecast for 2 year ahead for insample data. Please review my Rstan program which I successfully done for one year ahead. I need help for 2 year ahead.

``````dat=read.table("Loss Reserve Data_R.txt"); st=1;
N=dim(dat)[1] #number of rows
colnames(dat)=c("Year","Y","poli","develop","calender")
y=dat[,2]
y1=y[dat\$develop==1] #vector for 18 policy years
y0=0.7*y1

L=60000 #length of posterior sample

u=rep(0,L)

cal=as.numeric(dat[,5])
burn=5000
burn1=burn+1
iters=60000

dat\$poli<-factor(dat\$poli)
dat\$develop<-factor(dat\$develop)
dat\$calender<-factor(dat\$calender)

poli=as.integer(dat\$poli) # rows i effect (policy year)
develop=as.integer(dat\$develop) # columns j effect (development year)
calender=as.integer(dat\$calender) # columns n effect (calender year)

I=nlevels(dat\$poli);I1=I-1; #number of levels of policy (i)
J=nlevels(dat\$develop);J1=J-1 #number of levels of development (j)
K=nlevels(dat\$calender);K1=K-1 #number of calender year effect (n)

P=J1+7 #no of par in the model

#begin Rstan program

trialmodelA<- "

//Generalized Beta 2 Distribution

functions{
real beta2_log(real x, real a, real mu, real p, real q) {
return log(fabs(a))+(a*p-1)log(x)+lgamma(p+q)-(ap)log((mutgamma(p)*tgamma(q))/(tgamma(p+(1/a))*tgamma(q-(1/a))))-lgamma(p)-lgamma(q)-(p+q)log(1+pow((xtgamma(p+(1/a))tgamma(q-(1/a)))/(mutgamma(p)*tgamma(q)),a));
}
}
//data block
data {
int<lower=1> N; // number of observations
//int<lower=1> L; // number of posterior samples
real y[N]; // amounts of claims
real<lower=1> cal[N];
int<lower=1> I; //number of policy (i)
int<lower=1> J; //number of develop (j)
int<lower=1> K; //number of calender (n)
int<lower=1> I1; //dimension of alpha
int<lower=1> J1; //dimension of beta
int<lower=1> K1; //dimension of delta
int<lower=1, upper=I> poli[N]; //policy id
int<lower=1, upper=J> develop[N]; //develop id
int<lower=1, upper=K> calender[N]; //calender id
real y0[I];
}

parameters {
real <lower=0> mu0; //mean effect
real beta[J1];
real<lower=0> delta;
real<lower=0> gamma1;
real<lower=0, upper=0.9> gamma2;
real<lower=0> a;
real<lower=0> p;
real<lower=0> q;
}

transformed parameters {
real mu[N];
real beta18;

beta18=-sum(beta[1:J1]);

for (i in 1:N) {
if (develop[i]==1 && poli[i]!=18) {mu[i]=exp(mu0+beta[develop[i]]+deltalog(cal[i])+gamma1log(y0[poli[i]]));} //case A with y0
if (develop[i]==18) {mu[i]=exp(mu0+beta18+deltalog(cal[i])+gamma1log(y[i-1])+gamma2log(mu[i-1]));} //case B with beta18
if (poli[i]==18) {mu[i]=exp(mu0+beta[develop[i]]+deltalog(cal[i])+gamma1log(y0[poli[i]]));} //case C with alpha18
if (calender[i]==18 && develop[i]!=18 && poli[i]!=18) {mu[i]=exp(mu0+beta[develop[i]]+deltalog(cal[i])+gamma1log(y[i-1])+gamma2log(mu[i-1]));} //case D, delta18
if (develop[i]!=1 && develop[i]!=18) {mu[i]=exp(mu0+beta[develop[i]]+deltalog(cal[i])+gamma1log(y[i-1])+gamma2*log(mu[i-1]));} //else
}
}

model {
mu0~ normal(0,1);
beta[1:J1] ~ normal(0,1);
delta ~ normal(0,1);
gamma1 ~ gamma(0.1,1);
gamma2 ~ gamma(0.9,1);
a ~ gamma(0.1,1);
p ~ gamma(0.1,1);
q ~ gamma(0.1,1);

for (i in 1:N){
y[i]~beta2(a, mu[i], p, q);
}
}

generated quantities {
real loglik[N];
real dev;
real muf[19];
real muf1[20];
//real yf[19];
//real eyf[19];

for (i in 1:N){
loglik[i]=log(fabs(a))+(a*p-1)log(y[i])+lgamma(p+q)-(ap)*log((mu[i]*tgamma(p)*tgamma(q))/(tgamma(p+(1/a))*tgamma(q-(1/a))))-lgamma(p)-lgamma(q)-(p+q)*log(1+pow((y[i]*tgamma(p+(1/a))*tgamma(q-(1/a)))/(mu[i]*tgamma(p)*tgamma(q)),a));
}

dev = -2*sum(loglik[]);

//Forecast calendar year 19

muf[1]=exp(mu0+beta18+deltalog(19)+gamma1log(y[18])+gamma2log(mu[18])); //y1,19, note Y_1,18=Y_18 need beta19 but not replace beta19 by beta18
muf[2]=exp(mu0+beta18+deltalog(19)+gamma1log(y[35])+gamma2log(mu[35])); //y2,18, note Y_2,17=Y_17+18=35
muf[3]=exp(mu0+beta[17]+deltalog(19)+gamma1log(y[51])+gamma2log(mu[51])); //y3,17, note Y_3,16=Y_35+16=51
muf[4]=exp(mu0+beta[16]+deltalog(19)+gamma1log(y[66])+gamma2log(mu[66])); //y4,16, note Y_4,15=Y_51+15=66
muf[5]=exp(mu0+beta[15]+deltalog(19)+gamma1log(y[80])+gamma2log(mu[80])); //y5,15, note Y_5,14=Y_66+14=80
muf[6]=exp(mu0+beta[14]+deltalog(19)+gamma1log(y[93])+gamma2log(mu[93])); //y6,14, note Y_6,13=Y_80+13=93
muf[7]=exp(mu0+beta[13]+deltalog(19)+gamma1log(y[105])+gamma2log(mu[105])); //y7,13, note Y_7,12=Y_93+12=105
muf[8]=exp(mu0+beta[12]+deltalog(19)+gamma1log(y[116])+gamma2log(mu[116])); //y8,12, note Y_8,11=105+11=116
muf[9]=exp(mu0+beta[11]+deltalog(19)+gamma1log(y[126])+gamma2log(mu[126])); //y9,11, note Y_9,10=116+10=126
muf[10]=exp(mu0+beta[10]+deltalog(19)+gamma1log(y[135])+gamma2log(mu[135])); //y10,10, note Y_10,9=126+9=135
muf[11]=exp(mu0+beta[9]+deltalog(19)+gamma1log(y[143])+gamma2log(mu[143])); //y11,9, note Y_11,8=135+8=143
muf[12]=exp(mu0+beta[8]+deltalog(19)+gamma1log(y[150])+gamma2log(mu[150])); //y12,8, note Y_12,7=143+7=150
muf[13]=exp(mu0+beta[7]+deltalog(19)+gamma1log(y[156])+gamma2log(mu[156])); //y13,7, note Y_13,6=150+6=156
muf[14]=exp(mu0+beta[6]+deltalog(19)+gamma1log(y[161])+gamma2log(mu[161])); //y14,6, note Y_14,5=156+5=161
muf[15]=exp(mu0+beta[5]+deltalog(19)+gamma1log(y[165])+gamma2log(mu[165])); //y15,5, note Y_15,4=161+4=165
muf[16]=exp(mu0+beta[4]+deltalog(19)+gamma1log(y[168])+gamma2log(mu[168])); //y16,4, note Y_16,3=165+3=168
muf[17]=exp(mu0+beta[3]+deltalog(19)+gamma1log(y[170])+gamma2log(mu[170])); //y17,3, note Y_17,2=168+2=170
muf[18]=exp(mu0+beta[2]+deltalog(19)+gamma1log(y[171])+gamma2log(mu[171])); //y18,2, note Y_18,1=170+1=171
muf[19]=exp(mu0+beta[1]+deltalog(19)+gamma1log(y0[18])); //y19,1, note Y_19,0=Y0, need y0[19] but not replace y0[19] by y0[18]

//calendar year 19+1 ??

}
"
datt = list(N=N,y=y,y0=y0,cal=cal,I=I,J=J,K=K,I1=I1,J1=J1,K1=K1,poli=poli,develop=develop,calender=calender)

fitA = stan(model_code = trialmodelA,
model_name = "Trial modelA",
data = datt, iter = iters,
warmup=burn,
chains = 1)

print(fitA,digits = 4)
``````