Risk Calc

// Possible solutions implemented in Risk Calc for some of the challenge problems

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

// RAMAS RISK CALC CODE

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

// some ancillary functions

function many() {many_='#'; return $1.many_}

procedure data() {n_='#'; n_=$1.n_=_args-1; _for i_=1 _to n_ do $1.i_ = arg(i_+1); $1=string(n_)+' values: ' + string($2) + ',...' }

function sum() {n_='#'; n_=$1.n_; S_=0; _for i_=1 _to n_ _do S_=S_+$1.i_; _return S_;};

function mean() _if (1<_kind($1)) _return(_mean($1)) _else _return sum($1)/many($1);

function average() _return sum($1)/many($1);

function sd() _if (1<_kind($1)) _return(_sd($1)) _else {M_=mean($1); n_='#'; n_=$1.n_; S_=0; _for i_=1 _to n_ _do S_=S_+($1.i_-M_)^2; return(_sqrt(S_/(many($1)-1)))}

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

// PROBLEM 1.0

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

data(c, 12.1, 6.45, 73, 24.6, 15.2, 44.3, 19.0)

data(m, 58.22, 57.25, 53.23, 43.91, 42.69, 46.23, 43.08, 46.06, 44.6, 47.55, 53.91, 63.05, 43.62, 48.26, 62.83, 41.3, 39.05, 38.66, 55.3, 45.19, 68.7)

imean = 2.5

istd =0.5

///////////////////////////////////////////////////////////////////

// Method of matching moments, assuming independence

C = lognormal(mean(c), std(c))

M = normal(mean(m), sd(m))

I = normal(imean, istd)

H = C |*| I |/| M // vertical lines around operators indicate that we are assuming the operands are independent

5 < H // probability that 5 < H

///////////////////////////////////////////////////////////////////

// Confidence boxes, assuming independence

C = CBlognormal(c)

M = CBnormal(m)

I = normal(imean, istd)

H = C |*| I |/| M // vertical lines around operators indicate that we are assuming the operands are independent

5 < H // probability that 5 < H

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

// PROBLEM 1.1

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

C = env(0, CBlognormal(c))

M = CBnormal(m)

I = normal(imean, istd)

H = C |*| I |/| M // vertical lines around operators indicate that we are assuming the operands are independent

5 < H // probability that 5 < H

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

// PROBLEM 1.2

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

C = env(0, CBlognormal(c))

for i = 1 to many(m) do m.i = censor(m.i, [0,20])

M = CBnormal(m)

I = normal(imean, istd)

H = C |*| I |/| M // vertical lines around operators indicate that we are assuming the operands are independent

5 < H // probability that 5 < H

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

// PROBLEM 1.3

///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////

C = env(0, CBlognormal(c))

for i = 1 to many(m) do m.i = censor(m.i, [0,20])

M = CBnormal(m)

I = normal(imean, istd)

setdefault(copula, 'normal')

H = C * convolve(I, M, 0.4, '/') // intake and mass are assumed to be correlated with r=0.4

5 < H // probability that 5 < H