Plot various Normal's

[ac_1]

Hi Tom,

I would like to plot various theoretical Normal's
(i) individually
(ii) all superimposed
without having to specify the parameters in GRIDSERIES, i.e. have the plots automatically scaled on the x-axis. (The reason being I want to use in a procedure and pass in the mean's & sd's.)

I can do as follows e.g. plot 4 analytical PDF's as in (i) and (ii), but have had to specify @GRIDSERIES(from=-6,to=+6,n=1000,pts=ngrid) xsr

Code: Select all

* Various Normal's using PDF FORMULA
* - N(1.0,(0.1))
* - N(0.5,(1.1))
* - N(1.0,(3.0))
* - N(-1.2,(0.4))

@gridseries(from=-6,to=+6,n=1000,pts=ngrid) xsr
comp ac_pi = 3.141592653589793115998

comp mu1 = 1.0
comp sigma1 = sqrt(0.1)
set dn1 / = ( 1.0/(sigma1*sqrt(2.0*ac_pi)) * exp(-((xsr-mu1)^2.0)/(2.0*sigma1^2.0)) ); * or
set dn1 / = 1.0/sigma1*%density((xsr-mu1)/sigma1)

comp mu2 = 0.5
comp sigma2 = sqrt(1.1)
set dn2 / = ( 1.0/(sigma2*sqrt(2.0*ac_pi)) * exp(-((xsr-mu2)^2.0)/(2.0*sigma2^2.0)) ); * or
set dn2 / = 1.0/sigma2*%density((xsr-mu2)/sigma2)

comp mu3 = 1.0
comp sigma3 = sqrt(3.0)
set dn3 / = ( 1.0/(sigma3*sqrt(2.0*ac_pi)) * exp(-((xsr-mu3)^2.0)/(2.0*sigma3^2.0)) ); * or
set dn3 / = 1.0/sigma3*%density((xsr-mu3)/sigma3)

comp mu4 = -1.2
comp sigma4 = sqrt(0.4)
set dn4 / = ( 1.0/(sigma4*sqrt(2.0*ac_pi)) * exp(-((xsr-mu4)^2.0)/(2.0*sigma4^2.0)) ); * or
set dn4 / = 1.0/sigma4*%density((xsr-mu4)/sigma4)

spgraph(hea='',vfi=2,hfi=2)
scatter(key=upleft,style=line) 1
# xsr dn1 / 1
scatter(key=upleft,style=line) 1
# xsr dn2 / 2
scatter(key=upleft,style=line) 1
# xsr dn3 / 3
scatter(key=upleft,style=line) 1
# xsr dn4 / 4
spgraph(done)

scatter(key=upleft,style=line) 4
# xsr dn1
# xsr dn2
# xsr dn3
# xsr dn4

Alternatively, for each mean & sd pair i.e. a PDF, I can generate 1000000 random Normal variates, and using DENSITY and SCATTER, RATS automatically scales the PDF's on the plots; but the problem being as the PDF's are being 'simulated' the tails do not show all the way to LHS or RHS of the x-axis.

thanks,
Amarjit

[TomDoan]

You don't do anything to set the length of the workspace. Add the ALLOCATE instruction after @gridseries.

@gridseries(from=-6,to=+6,n=1000,pts=ngrid) xsr
all ngrid

"pi" is %PI. (Though obviously the %DENSITY calculation is what you should use).

[ac_1]

Thanks.

The analytical normal PDF's look almost identical as simulating 1000000 random normal variates! But the scaling on the x-axis is fixed with former, which is better for comparative purposes, i.e.

@gridseries(from=(minmean-(3.0*maxsigma)),to=(maxmean+(3.0*maxsigma)),n=1000,pts=ngrid) xsr
all ngrid

I can do the same for the log-normal distribution .

But cannot for the non-central chi-square (which is the back-transformed distribution from a SQRT transformation).

From https://estima.com/ratshelp/index.html? ... ntral.html, the PDF is a complicated formula with parameters r>0 (degrees of freedom), and nc>0 (non-centrality parameter). For the normal plots example

mean sigma
1.0, sqrt(0.1)
0.5, sqrt(1.1)
1.0, sqrt(3.0)
-1.2, sqrt(0.4)

how to plot the distributions using %CHISQRNCDENSITY(x,r,nc) i.e. what are x, r and nc, for the 4 means's and sigma's?

[TomDoan]

Thanks.

The analytical normal PDF's look almost identical as simulating 1000000 random normal variates! But the scaling on the x-axis is fixed with former, which is better for comparative purposes, i.e.

@gridseries(from=(minmean-(3.0*maxsigma)),to=(maxmean+(3.0*maxsigma)),n=1000,pts=ngrid) xsr
all ngrid

I can do the same for the log-normal distribution .

But cannot for the non-central chi-square (which is the back-transformed distribution from a SQRT transformation).

From https://estima.com/ratshelp/index.html? ... ntral.html, the PDF is a complicated formula with parameters r>0 (degrees of freedom), and nc>0 (non-centrality parameter). For the normal plots example

mean sigma
1.0, sqrt(0.1)
0.5, sqrt(1.1)
1.0, sqrt(3.0)
-1.2, sqrt(0.4)

how to plot the distributions using %CHISQRNCDENSITY(x,r,nc) i.e. what are x, r and nc, for the 4 means's and sigma's?

There are formulas in the help link for mean and variance (not standard deviation) based upon the shape parameters. Backsolve for the underlying parameters.

[ac_1]

Thanks.

For some means and sd's (mu must be integer and delta must be non-negative) there will be no solution.

E.g. for mean=1, sd=sqrt(0.1)
mu + delta = 1
2*mu + 4*delta = 0.1

But e.g. mean=16, sd=sqrt(70) does have a solution
mu + delta = 16
2*mu + 4*delta = 70
mu=-3, delta=19

Hence it's best I use the 1000000 simulation method to plot various non-central chi-squared's.


Also, after transforming a series => to back-transform, are these the correct distributions?

NONE [x^1] => [x^1] = normal
SQRT [x^0.5] => [x^2] = x^1 = non-central chi-squared
LOG [log(x)] => [exp(x)] = x^1 = log-normal --- why is it named the log-normal and not exp-normal?
INVERSE [x^-1] => [x^2] = x^1 = non-central chi-squared
CUBE-ROOT [x^(1/3)] => [x^3] = x^1 = ID?

And, what are the probability functions PDF, CDF, Inverse, Random draws, etc, for the Box-Cox distributions (Power-Normal distributions)?
https://en.wikipedia.org/wiki/Box%E2%80 ... tribution.

[TomDoan]

I think you're forgetting that in your application, your degrees of freedom is 1 (you're just squaring a single random variable), but it's a scaled non-central chi-squared. So you need to peg nu at 1, and compute the delta and scale factor which gives you the desired mean and standard deviation (or variance). And yes, there are some mean/variance combinations which are incompatible with that---I'm not sure what you think "simulating" will do, as if it's impossible for a scaled chi-squared(1) to produce a particular ratio of mean to standard deviation, there is no way to simulate one which has that combination.

[ac_1]

Sorry, typo above, nu not mu

I think you're forgetting that in your application, your degrees of freedom is 1 (you're just squaring a single random variable), but it's a scaled non-central chi-squared. So you need to peg nu at 1, and compute the delta and scale factor which gives you the desired mean and standard deviation (or variance). And yes, there are some mean/variance combinations which are incompatible with that

I am unclear how to do this in RATS, and automate the backsolve for the underlying parameters for a series of mean's and std's. E.g let's say the results from forecasting the SQRT of a series are:

Code: Select all

mean            std
128.480485268   0.459514910574
128.659193181   0.630766425674
128.849570264   0.755998905248
129.020870238   0.859639425665
129.201289182   1.149002767719
129.380784764   1.233648748511

the mean's are then SQUARED into the original scale.

I'm not sure what you think "simulating" will do, as if it's impossible for a scaled chi-squared(1) to produce a particular ratio of mean to standard deviation, there is no way to simulate one which has that combination.

I can plot a theoretical a non-central chi squared, e.g. as

Code: Select all

@gridseries(from=0,to=50,n=1000,pts=ngrid) xsr1
all ngrid
comp nu1 = 4
comp delta1 = 1
set dncchisq1 / = %CHISQRNCDENSITY(xsr1,nu1,delta1)
scatter(key=upleft,style=line) 1
# xsr1 dncchisq1

or I can simulate various non-central chi squared's, given a series of mean's and std's in the SQRT scale, e.g. as above

Code: Select all

comp ndraws = 1000000
do h = 1, nSteps
	set dtn(h) 1 ndraws = (mean(startl+h-1)+std(startl+h-1)*%ran(1.0))**2.0
end do h

and plot with DENSITY and SCATTER, back into the original scale.

The aim ideally being in the original scale, plot 6 theoretical PDF's in this example (knowing values for mean, std, and backsolving for delta, scale, pegging nu=1), in preference to simulate 1000000 times (where I only need the mean's and std's).

[TomDoan]

It's relatively simple algebra. If you multiply X by c, the mean changes by a factor of c, the variance by a factor of c^2. nu is fixed at 1. You have to solve for c and the non-centrality parameter delta. You can eliminate c by looking at mu^2/sigma^2 and solving that for delta. Given that, solve for c.

[ac_1]

From https://estima.com/ratshelp/index.html? ... ntral.html

nu + delta = mean
2*nu + 4*delta = sigma^2

It's relatively simple algebra. If you multiply X by c, the mean changes by a factor of c, the variance by a factor of c^2. nu is fixed at 1. You have to solve for c and the non-centrality parameter delta. You can eliminate c by looking at mu^2/sigma^2 and solving that for delta. Given that, solve for c.

nu + delta = mean*c
2*nu + 4*delta = sigma^2*c^2

Fix the degrees of freedom nu = 1

1 + delta = mean*c
2*1 + 4*delta = sigma^2*c^2

delta = mean*c - 1
4*delta = sigma^2*c^2 - 2

Put c = mean^2/sigma^2

delta = mean*(mean^2/sigma^2) - 1
4*delta = sigma^2*(mean^2/sigma^2)^2 - 2

delta = mean*(mean^2/sigma^2) - 1
4*delta = sigma^2*(mean^4/sigma^4) - 2

delta = mean*(mean^2/sigma^2) - 1
4*delta = (mean^4/sigma^2) - 2

delta = (mean^3/sigma^2) - 1
4*delta = (mean^4/sigma^2) - 2

I know the mean & sigma, so

Put A = (mean^3/sigma^2)
Put B = (mean^4/sigma^2)

Therefore
delta = A - 1
4*delta = B - 2

Now, add the 2 equations, solving for delta

5*delta = A + B - 1 - 2

5*delta = A + B - 3

Therefore delta = (A + B - 3)/5

Hence, subsitutute in A & B, the underlying parameters are:
delta = ((mean^3/sigma^2) + (mean^4/sigma^2) - 3)/5
nu = 1

and c = mean^2/sigma^2

Correct

[TomDoan]

No. c is one of the unknowns; you can't just say it's equal to something.

nu + delta = mean*c
2*nu + 4*delta = sigma^2*c^2

As I described, square the first and divide to eliminate the c.

[ac_1]

Sorry I don't understand exactly what you are saying.

square the first

Does that mean square each side of

nu + delta = mean*c

if so,

(nu + delta)^2 = (mean*c)^2

nu^2 + 2*nu*delta + delta^2 = mean^2*c^2

and divide to eliminate the c.

nu^2 + 2*nu*delta + delta^2 = mean^2*c^2
2*nu + 4*delta = sigma^2*c^2

dividing is

0.5*nu + 0.5*nu + delta^2 = mean^2/sigma^2

nu + delta^2 = mean^2/sigma^2

solving for delta

delta = sqrt(mean^2/sigma^2 - nu)

peg nu = 1

delta = sqrt(mean^2/sigma^2 - 1)


But that's not right because if I try to plot the 1st PDF, mean and sigma from forecasting a SQRT series:

Code: Select all

@gridseries(from=0,to=500,n=1000,pts=ngrid) xsr1
all ngrid
comp nu1 = 1
comp mean1 = 128.480485268   
comp sigma1 = 0.459514910574
comp delta1 = sqrt((mean1^2/sigma1^2)-nu1)
set dncchisq1 / = %CHISQRNCDENSITY(xsr1,nu1,delta1)
scatter(key=upleft,style=line) 1
# xsr1 dncchisq1

the plot should be in SQ’s i.e. the original scale.

[TomDoan]

If A=B, then A^2=B^2. If C=D then A^2/C=B^2/D. In this case, B^2 has a c^2 factor, as does D. So B^2/D eliminates the c^2 leaving you with a (quadratic) equation for delta.

[ac_1]

If A=B, then A^2=B^2. If C=D then A^2/C=B^2/D. In this case, B^2 has a c^2 factor, as does D. So B^2/D eliminates the c^2 leaving you with a (quadratic) equation for delta.

What is the quadratic equation for delta in terms of mean, sigma, nu=1 (square of the standard normal), c has been eliminated, so that I may plot the theoretical PDF's?

[TomDoan]

A^2 is quadratic in delta, C is linear in delta. B^2/D is known since the unknown c^2 has cancelled. Multiply the equation up by C and convert to a quadratic equation.

At any rate, I assume what you are trying to do is to take an estimated density which is Normal in sqrt(x) and convert to a density for X. You don't need to force that into a non-central chi-squared since you can just do the density of the transformed variable. If you aren't familiar with this, there's a primer at

https://www.cl.cam.ac.uk/teaching/0708/ ... prob11.pdf

It's basically the same thing that is done to derive the density of the non-central chi-squared and the log normal, but more general for the chi-squared (since it doesn't assume a unit variance).

This is an example where the sqrt(x) is assumed N(mu,sigma^2). The fx is computing the density of X using the fact that the density of sqrt(X) is known---you use the known density inserting sqrt(X), and multiply by the derivative of sqrt(X) with respect to X.

*
* Mean and std dev of sqrt(x)
*
compute mu=3.0,sigma=0.5
@gridseries(from=0.0,to=20.0,size=0.1) x
set fx = exp(%logdensity(sigma^2,sqrt(x)-mu))/(2*sqrt(x))
scatter
# x fx

I'm not sure what you expect to find from graphing those. They will have fairly predictable behaviors. A sqrt transformation will invert to a slightly heavier right tail. (A log transformation will be heavier still). That's basically by construction.

[ac_1]

Thanks!

I'm not sure what you expect to find from graphing those. They will have fairly predictable behaviors. A sqrt transformation will invert to a slightly heavier right tail.
(A log transformation will be heavier still). That's basically by construction.

Yes.

Plotted as a "side-view" of fan charts, as forecast PDF's at EACH time-step ahead. A fan chart as a time plot IMHO can be improved upon as it's not the whole picture. A side-view provides a useful extra-dimension, especially with superimposed PDF's as comparisons between the heights of the densities can easily be seen e.g. typically in multi-step ahead forecasts, where the PDF's at near horizons will be tall and thin, and wider and shorter at longer horizons.

In an ideal world I would like a 3D fan chart in RATS, as p.35 here https://www.bankofengland.co.uk/-/media ... -chart.pdf

Interestingly: can a probability be calculated as the difference between the shapes of the heights? i.e. the area under a PDF=1, if I have 2 PDF's with different shapes can the "gaps" be taken into account as the difference in probabilities, or is that the log-predictive density where the maximum height or density is calculated and averaged across ALL forecasts, but is only applicable with Gaussian PDF's? The idea being the forecast at time-step 1 has a higher probability rather than density - a good explanation here: https://math.stackexchange.com/question ... 29%20d%20x., than the forecast at time-step 2, time-step 3, etc etc.


What about?

Also, after transforming a series => to back-transform, are these the correct distributions?

NONE [x^1] => [x^1] = normal
SQRT [x^0.5] => [x^2] = x^1 = non-central chi-squared
LOG [log(x)] => [exp(x)] = x^1 = log-normal --- why is it named the log-normal and not exp-normal?
INVERSE [x^-1] => [x^2] = x^1 = non-central chi-squared
CUBE-ROOT [x^(1/3)] => [x^3] = x^1 = ID?

And, what are the probability functions PDF, CDF, Inverse, Random draws, etc, for the Box-Cox distributions (Power-Normal distributions)?
https://en.wikipedia.org/wiki/Box%E2%80 ... tribution.