inst/doc/intro.R

In Jinshuyue98/SC190: Final homework for 'Statistical Computing' course

## ------------------------------------------------------------------------
Equity_after_shock <- function(Epreall,Aall,Lall,Aball,Lball,p,k,a,b,beta,R){
 # Epreall the equity before the shock
 # Aall total assets
 # Lall total liabilities
 # Aball inter-bank assets
 # Lball inter-bank liabilities
 # p probability of inter-bank liabilities
 # k parameter modelling the capital cushion
 # a,b determining the shape of the cdf of the Beta distribution
 # beta parameter modelling the actual exogenous recovery rate
 # R parameter modelling the perceived exogenous recovery rate 

  col <- ncol(Aall)
  row <- nrow(Aall)
  Eeall<-matrix(data=NA, nrow = row, ncol = col)
  for (h in 1:col){
    Epre<-Epreall[,h]
    Ee<-Epre
    A<-Aall[,h]
    L<-Lall[,h]
    Ab<-Aball[,h]
    Lb<-Lball[,h]
    Ae<-A-Ab
    x<-0.05*Ae
    Absum<-sum(Ab[])
    lamda<-p*col*(col-1)/Absum
    nsamples<-5
    s<-sample_ERE(Lb, Ab, p, lamda, nsamples, thin = 5000, burnin = 10000)
    ssum<-s[[1]]
    for (i in 2:nsamples) {ssum<-ssum+s[[i]]}
    Lij<-ssum/nsamples
    g<-1
    for (g in 1:1000){
      ab<-numeric(row)
      i<-1
      j<-1
      for (i in 1:row) {
        for (j in 1:row) {
          y<-(Ee[j]+L[j])/L[j]
          if (y>=1+k) ab[i]<-ab[i]+Lij[j,i]
          if (y>=1&&y<1+k) {
            F<-pbeta((1+k-y)/k,a,b)
            ab[i]<-ab[i]+Lij[j,i]*(1-(1-R)*F)}
          if (y>=0&&y<1) ab[i]<-ab[i]+Lij[j,i]*beta*y
          if (y<0) ab[i]<-ab[i]
        }
        Ee[i]<-Ae[i]-x[i]+ab[i]-L[i]
      }
    }
    Eeall[,h]<-Ee
  }
  return(Eeall)
   # a equity matrix of the banking network after the shock
}

## ------------------------------------------------------------------------
A2018<-c(2769954000,2322269300,2126727500,2260947100,953117100,674572900)
L2018<-c(2535465700,2123109900,1954187800,2093468400,882586300,620212400)
Epre2018<-c(234488300,199159400,172539700,167478700,70530800,54360500)
Ab2018<-c(57780300,34972700,78176100,55201300,56477800,31341100)
Lb2018<-c(70461833,60893785,47421310,47173804,58444388,29554180)
A2017<-c(2608704300,2212438300,1946742400,2105338200,903825400,629763800)
L2017<-c(2394598700,2032855600,1789074500,1962398500,836198300,581424600)
Epre2017<-c(214105600,179582700,157667900,142939700,67627100,48339200)
Ab2017<-c(47753700,32523300,48655900,50526900,57076600,15462800)
Lb2017<-c(58630225,45721989,28804785,33377591,52960249,32504361)
A2016<-c(2413726500,2096370500,1814888900,1957006100,840316600,594231100)
L2016<-c(2215610200,1937405100,1666179700,1824847000,777075900,553894900)
Epre2016<-c(198116300,158965400,148709200,132159100,63240700,40336200)
Ab2016<-c(52741500,26067000,48392900,58094900,46779100,20025100)
Lb2016<-c(65863017,42478615,24550718,39775516,46656201,32776433)
A2015<-c(2220978000,1834948900,1681559700,1779139300,715536200,547497800)
L2015<-c(2040926100,1690440600,1545799200,1657950800,661727000,511322000)
Epre2015<-c(180051900,144508300,135760500,121188500,53809200,36175800)
Ab2015<-c(47223400,31077900,35021800,50425200,35681200,18569300)
Lb2015<-c(57837485,38959978,32024949,38239057,29287800,21649531)
A2014<-c(2060995300,1674413000,1525138200,1597415200,626829900,473182900)
L2014<-c(1907264900,1549176700,1406795400,1494153300,579469400,441676900)
Epre2014<-c(153730400,125236300,118342800,103261900,47360500,31506000)
Ab2014<-c(47850300,24852500,29911100,40706200,17231800,12408500)
Lb2014<-c(55175483,25823315,24020165,28696639,27174942,12069856)
A2013<-c(1891775200,1536321000,1387429900,1456210200,596093700,401639900)
L2013<-c(1763928900,1428888100,1291282200,1371756500,553945300,375044300)
Epre2013<-c(127846300,107432900,96147700,84453700,42148400,26595600)
Ab2013<-c(41161800,15206500,40200100,30865500,19705700,14804700)
Lb2013<-c(48795965,18918097,32506094,21156230,25385098,15182816)
A2012<-c(1754221700,1397282800,1268061500,1324434200,527337900,340821900)
L2012<-c(1641375800,1302321900,1181907300,1249298800,489193200,320771200)
Epre2012<-c(112845900,94960900,86154200,75135400,38144700,20050700)
Ab2012<-c(22451300,12965300,34925100,22338000,15378700,10342000)
Lb2012<-c(27821819,13162776,26657313,16387905,22350645,12019942)
A2011<-c(1547686800,1228183400,1183006600,1167757700,461117700,279497100)
L2011<-c(1451904500,1146517300,1107417200,1102778900,433838900,262996100)
Epre2011<-c(95782300,81666100,75589400,64978800,27278800,16501000)
Ab2011<-c(16051600,10904000,45466800,21268000,14872600,13138100)
Lb2011<-c(33765670,10641493,31251723,14727772,22192429,9122013)
A2010<-c(1345862200,1081031700,1045986500,1033740600,395159300,240250700)
L2010<-c(1263696500,1010941200,978371500,979517000,372793600,226850100)
Epre2010<-c(82165700,70090500,67615000,54223600,22365700,13400600)
Ab2010<-c(6491800,6396200,13698400,9537500,7961000,5891700)
Lb2010<-c(11457752,6044018,14186825,5171232,8960508,4156265)
A2009<-c(1178505300,962335500,875194300,888258800,330913700,206794100)
L2009<-c(1110611900,906433500,820654900,853966300,314471200,197515800)
Epre2009<-c(67893400,55902000,54539400,34292500,16442500,9278300)
Ab2009<-c(7790600,2221700,13893300,4943500,4329700,6239700)
Lb2009<-c(7526027,4062247,11167556,2803930,9135686,4723054)
A2008<-c(975765400,755545200,695569400,701435100,267825500,157179700)
L2008<-c(915051600,708789000,646179300,672381000,253261300,149201600)
Epre2008<-c(60713800,46756200,49390100,29054100,14564200,7978100)
Ab2008<-c(12679200,1683600,38574800,4447900,8953900,8183600)
Lb2008<-c(9565698,7686508,35243319,6085835,9194103,6747537)
Aall<-cbind(A2008,A2009,A2010,A2011,A2012,A2013,A2014,A2015,A2016,A2017,A2018)
Lall<-cbind(L2008,L2009,L2010,L2011,L2012,L2013,L2014,L2015,L2016,L2017,L2018)
Epreall<-cbind(Epre2008,Epre2009,Epre2010,Epre2011,Epre2012,Epre2013,Epre2014,Epre2015,Epre2016,Epre2017,Epre2018)
Aball<-cbind(Ab2008,Ab2009,Ab2010,Ab2011,Ab2012,Ab2013,Ab2014,Ab2015,Ab2016,Ab2017,Ab2018)
Lball<-cbind(Lb2008,Lb2009,Lb2010,Lb2011,Lb2012,Lb2013,Lb2014,Lb2015,Lb2016,Lb2017,Lb2018)

library("systemicrisk")
p<-0.5
k<-0.025
a<-b<-1
beta<-0.2
R<-0.5
Equity_after_shock(Epreall,Aall,Lall,Aball,Lball,p,k,a,b,beta,R)

## ------------------------------------------------------------------------
CI_BCa_Bootstrap <-function(x, conf) {
  # x the sample
  # conf confidence level

  x <- as.matrix(x)
  n <- nrow(x) 
  N <- 1:n
  alpha <- (1 + c(-conf, conf))/2
  zalpha <- qnorm(alpha)
  cov.hat<-cov(x)
  lamda.hat<-eigen(cov.hat)$values
  lamda.sum<-sum(lamda.hat)
  th0<-lamda.hat[1]/lamda.sum
  B<-2000
  th<-numeric(B)
  for(b in 1:B){
    i <- sample(1:n, size = n, replace = TRUE)
    scor.boot<-scor[i,]
    cov.boot<-cov(scor.boot)
    lamda.boot<-eigen(cov.boot)$values
    th[b]<-lamda.boot[1]/sum(lamda.boot)
  }
  z0 <- qnorm(sum(th < th0) / length(th))
  scor.jack<-x[-1,]
  theta.jack<-numeric(n)
  for(i in 1:n){
    scor.jack<-scor[-i,]
    cov.jack<-cov(scor.jack)
    lamda.jack<-eigen(cov.jack)$values
    theta.jack[i]<-lamda.jack[1]/sum(lamda.jack)
  }
  L <- mean(theta.jack) - theta.jack
  a <- sum(L^3)/(6 * sum(L^2)^1.5)
  adj.alpha <- pnorm(z0 + (z0+zalpha)/(1-a*(z0+zalpha)))
  limits <- quantile(th, adj.alpha, type=6)
  return(list("est"=th0, "BCa"=limits))
    # a list including the sample estimate and BCa confidence interval
}

## ------------------------------------------------------------------------
library("bootstrap")
data(scor,package="bootstrap")
CI_BCa_Bootstrap(scor,conf=0.95)

## ----eval=FALSE----------------------------------------------------------
#  NumericVector Metropolis(double sigma, double x0, int N){
#    #x0: the initial point
#    #sigma: the standard deviation in the normal distribution
#    #N: the length of the chain
#    NumericVector x(N);
#    x[0]=x0;
#    for (int i=1;i<N;i++) {
#      double e=runif(1)[0];
#      double z=rnorm(1,x[i-1],sigma)[0];
#      if (e<=exp(abs(x[i-1])-abs(z))) x[i]=z;
#      else {
#        x[i]=x[i-1];
#      }
#    }
#    return x;
#  }

## ----eval=TRUE-----------------------------------------------------------
library(SC19083)
N<-100
sigma<-1.0
x0<-0.0
Metropolis(sigma,x0,N)