Plot various Normal's

[TomDoan]

I find most 3-D graphs used in economics and finance to be busy and not particularly helpful. (Basically, the viewer is forced to try to visualize a cut into two dimensions anyway).

The question about percentiles is why I asked what you think the density graphs are giving you. People cannot visually integrate which is why you can't easily look at a density graph and tell much of anything about it---you can pick out the mode and (if graphed in a pair) tell which has fatter tails. That's it. And the mode doesn't really say much, particularly for skewed distributions. Look at densities of chi-squared for various degrees of freedom: the mode for low DF is zero, while for slightly higher DF, the density at zero is 0.

Regardless of which monotonic transformation you use, if g(X)~N(mu,sigma^2), then the percentiles of X can be obtained by inverting the corresponding percentiles of the Normal. So you don't need the density function---you just need to know the inverse function for g. And you don't need a name for a distribution since you can always just do the transformed density of the Normal. If you need it. Which you probably don't.

[ac_1]

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

The example above has sqrt(x) assumed N(mu,sigma^2), i.e. mu is the mean.


Just to be clear:

mu=3.0 and sigma=0.5 are forecasts from a SQRT transformation.

To back-transform the point forecast I simply SQUARE mu, that will be the MEDIAN: mu^2.

The back-transformed MEAN can be calculated from a Taylor Series expansion: mu^2 + sigma^2.


For comparison here's the 1000000 simulation method.

comp ndraws = 1000000
set mu1 1 ndraws = 3.0
set sigma1 1 ndraws = 0.5
set dtn 1 ndraws = (mu1+sigma1*%ran(1.0))**2.0
density(smoothing=2.0) dtn 1 ndraws xsr mu1
scatter(key=upright,footer="Conditional PDFs",style=line,hmax=20) 1
# xsr mu1 1 ndraws 1


The plots are almost identical. The simulation: is automatically scaled (although I have hmax=20 to resemble the first); the theoretical: the limits have to be determined in @gridseries.

I think they are both correct?

[TomDoan]

Other than the fact that the mean of the transformed variable isn't obtained by a Taylor expansion, but simply because EX^2=(EX)^2+var(X) in general. For other transformations and other distributions, there is no relationship between the moments you get from a Taylor expansion and the actual moments.

[ac_1]

Thanks - on the topic of comparing distributions and referring to https://www.estima.com/textbooks/johnst ... n4p351.rpf, there's a typo

compute b(draw)=%beta(1),tstat(draws)=%tstats(1)

should be

compute b(draw)=%beta(1),tstat(draw)=%tstats(1)

If seed 1, then from stats(fractiles) b
Sample Mean 0.005308
Skewness 92.551084
Median 0.004980

As to be expected the Mean is to the right of the Median for positively skewed data. The MODE (the maximum height) looks to be at 0.0051.

If seed 500, then from stats(fractiles) b
Sample Mean 0.004986
Skewness -14.950009
Median 0.004981

But not as expected for negatively skewed data, the Mean is to the right of the Median, as compared with
https://www.bing.com/images/search?view ... ajaxserp=0

and the MODE appears to be at 0.00505, a bit less.

So I cannot just look at Skewness to ascertain the positions of the Mean and Median of any distribution e.g. let's say for a specific forecast density at step h via bootsrapping. Why?

Also,

I tried to calculate the MODE of a continuous distribution (sample) which would be the max of the density, example here https://en.wikipedia.org/wiki/Mode_(statistics), but don't know how in RATS?

And further, is there a measure similar to LPDS (only for a Gaussian PDF and using the Mean) but using the MODE to assess the accuracy of a fan (forecast MODE, forecast std and actual), which can be averaged across EACH time-step for ANY forecast distribution, empirical and theoretical?

[TomDoan]

Re the johnston example, the version of that which is distributed with RATS has been corrected for many years. We'll fix the one on the web site, but it's unlikely very many people have ever looked for that. (It's an old textbook).

That is a *very* misleading comparison of distributions. (I'm not sure where you got that, but the wikipedia description is more accurate). The definition of skew is *NOT* that the median is to a particular side of the mean; it's that the central third moment is non-zero. 2SLS has very fat tails (see the size of the kurtosis) which means that there are a few particularly large values. If, in a simulation, those really big values are mainly negative, the skewness will be negative; if they are positive, the skewness will be positive. Even though the median is scarcely affected. (In reality, the distribution is symmetric, so there is no "skew"; it's just that the 3rd and 4th moments don't exist, so the sample values for those never really settle down).

If you have an actual parametric distribution, you can do FIND MAX to find the mode. If you have an empirical distribution, you can do a DENSITY instruction and approximate the mode by doing an EXTREMUM to find the maximum value.

I doubt very seriously that there is anything that would analyze the modes of forecasts. There is no sensible loss function which would ever lead you to use the mode of a distribution as a central estimate. (Median minimizes absolute deviations, mean minimizes squared deviations). The only real advantage of the mode is that you can identify it from a graph.

[ac_1]

Re the johnston example, the version of that which is distributed with RATS has been corrected for many years. We'll fix the one on the web site, but it's unlikely very many people have ever looked for that. (It's an old textbook).

https://www.estima.com/cgi-bin/bookbrowser.cgi I searched Topic Distribution (empirical).

That is a *very* misleading comparison of distributions. (I'm not sure where you got that, but the wikipedia description is more accurate). The definition of skew is *NOT* that the median is to a particular side of the mean; it's that the central third moment is non-zero. 2SLS has very fat tails (see the size of the kurtosis) which means that there are a few particularly large values. If, in a simulation, those really big values are mainly negative, the skewness will be negative; if they are positive, the skewness will be positive. Even though the median is scarcely affected. (In reality, the distribution is symmetric, so there is no "skew"; it's just that the 3rd and 4th moments don't exist, so the sample values for those never really settle down).

The comparson is from University days.

If you have an actual parametric distribution, you can do FIND MAX to find the mode. If you have an empirical distribution, you can do a DENSITY instruction and approximate the mode by doing an EXTREMUM to find the maximum value.

I doubt very seriously that there is anything that would analyze the modes of forecasts. There is no sensible loss function which would ever lead you to use the mode of a distribution as a central estimate. (Median minimizes absolute deviations, mean minimizes squared deviations). The only real advantage of the mode is that you can identify it from a graph.

Code: Select all

* parametric distribution

@gridseries(from=0,to=+8,n=1000,pts=ngrid) xsr
all ngrid

comp mu = 4.0
comp sigma = 0.75
set dnorm / = 1.0/sigma*%density((xsr-mu)/sigma)
*
prin / xsr dnorm

scatter(style=lines) 1
# xsr dnorm

* convert series to array
make adnorm
# dnorm

* max y-value
comp maxarray = %maxvalue(adnorm)
disp maxarray

I can calculate the max y-value i.e. height 0.53192, but not the x-value i.e. MODE.

(For a standard normal mu = 0, sigma = 1.0, the MODE on the x-axis is obviously 0, y-axis 0.39894, also via prin / xsr dnorm)

To use FIND how do I convert the adnorm array RECTANGULAR[REAL] to a REAL i.e. calculating the max height with the associated x-value MODE?

Code: Select all

dec real dens
find maximum dens
   comp dens = adnorm
end find

Code: Select all

* empirical distribution

seed 71271
comp ndraws = 1000
set mu 1 ndraws = 5.0
set sigma 1 ndraws = 0.5
set dtn 1 ndraws = (mu+sigma*%ran(1.0))

stats(fractiles) dtn

density(smoothing=2.0) dtn 1 ndraws xs mu
scatter(footer="",style=line,hgrid=||%mean,%median||) 1
# xs mu 1 ndraws 1

extremum dtn

From eye the MODE on the x-axis is approx. 4.9, y-axis (max density, dtn) approx. 0.81

extremum dtn doesn't calculate the maximum height at the x-value?

I doubt very seriously that there is anything that would analyze the modes of forecasts. There is no sensible loss function which would ever lead you to use the mode of a distribution as a central estimate. (Median minimizes absolute deviations, mean minimizes squared deviations). The only real advantage of the mode is that you can identify it from a graph.

Similar to LPDS (Gaussian distribution only and for the mean), is there an equivalent loss measure for a LOG t/f (log-normal distribution) or for a SQRT t/f (non-central chi-squared distribution)?

[TomDoan]

FIND can be used for parametric densities (where you have a function you can evaluate at an arbitrary value). EXTREMUM is used for data series, where you are just trying to find the maximum observed value. EXTREMUM defines %MAXENT which is the entry at which the maximum occurs. Just evaluate the "x" series at that point.

The LPDS has nothing to do with Gaussianity---it can be used with whatever density you want. Again note (as I did when you first asked about it), that it measures how well the predictions match an assumed distribution, not how good the predictions are. (Worse predictions with a better variance estimate can lead to a better LPDS).

[ac_1]

EXTREMUM is used for data series, where you are just trying to find the maximum observed value. EXTREMUM defines %MAXENT which is the entry at which the maximum occurs. Just evaluate the "x" series at that point.

Here's my steps

Code: Select all

seed 71271
comp ndraws = 1000
set mu 1 ndraws = 5.0
set sigma 1 ndraws = 0.5
set dtn 1 ndraws = (mu+sigma*%ran(1.0))

View -> Series Window
MU 1000
SIGMA 1000
DTN 1000

Code: Select all

density(smoothing=2.0) dtn 1 ndraws xs mu

View -> Series Window
MU 100
SIGMA 1000
DTN 1000
XS 100

I've then put maxgrid=ndraws

Code: Select all

density(smoothing=2.0,maxgrid=ndraws) dtn 1 ndraws xs mu

View -> Series Window
MU 1000
SIGMA 1000
DTN 1000
XS 1000

And

Code: Select all

extremum mu

Extreme Values of Series MU
Minimum Value is 0.00296952612 at Entry 1000
Maximum Value is 0.80475188010 at Entry 423

So, the max density is 0.80475188010 at Entry 423

Now,

Code: Select all

order(index=ix) dtn
set odtn = dtn(ix)
prin / odtn

odtn entry 423 is 4.89470012320

So, the MODE is 4.89470012320 with density 0.80475188010.

If I plot

Code: Select all

stats(fractiles) dtn
scatter(footer="",style=line,hgrid=||%mean,%median,4.89470012320||) 1
# xs mu 1 ndraws 1

The MODE is not at the peak.

The LPDS has nothing to do with Gaussianity---it can be used with whatever density you want. Again note (as I did when you first asked about it), that it measures how well the predictions match an assumed distribution, not how good the predictions are. (Worse predictions with a better variance estimate can lead to a better LPDS).

These are the following calculations (in a procedure with if statements), I think they are correct. Please check.
Are they beneficial in a backtest summary/loss measures table?

Code: Select all

* NORMAL
set scores startl endl = ( 1.0/(fsigma*sqrt(2.0*%pi)) * exp(-((actual-fcast)^2.0)/(2.0*fsigma^2.0)) )
set lscores startl endl = log(scores)
sstats(mean) startl endl lscores>>LPDS
comp %%LPDS = LPDS

* LOG-NORMAL
set scores startl endl = ( 1.0/(actual*fsigma*sqrt(2.0*%pi)) * exp(-((log(actual)-fcast)^2.0)/(2.0*fsigma^2.0)) )
set lscores startl endl = log(scores)
sstats(mean) startl endl lscores>>LPDS
comp %%LPDS = LPDS

* NON-CENTRAL CHI-SQUARED
set scores startl endl = %logdensity(fsigma**2.0,actual-fcast)-log(2)-0.5*log(actual)
set lscores startl endl = log(scores)
sstats(mean) startl endl lscores>>LPDS
comp %%LPDS = LPDS

[TomDoan]

1. I don't know what "mu1" is. And why are you even using "mu" as the target series for the estimated density? It has nothing to do with the mean.

2. What's wrong with 100 ordinates for the density? There is no justification for those being the same as NDRAWS---you're not making anything more accurate by doing that. (The estimated density at each X(i) is computed across the full sample regardless of how many X's you choose to use).

3. I have no idea what you are trying to do with the ORDER instruction. You found the maximum value at 483. The "x" at that is just xs(483) (or more generally xs(%maxent)).

4. The mean and variance for the log normal and square root normal are not the same as the mean and variance for the normal itself; they are the mean and variance of the transformed data (or in your case the prediction and predictive variance of the transformed data). And it would be sqrt(actual) in the %logdensity for the latter.

[ac_1]

Re 1. 2. 3. This is the correct figure

Code: Select all

seed 71271
comp ndraws = 1000
set mu 1 ndraws = 5.0
set sigma 1 ndraws = 0.5
set dtn 1 ndraws = (mu+sigma*%ran(1.0))

density(smoothing=2.0) dtn 1 ndraws xs mu

prin / mu

extremum mu

comp mode = xs(%maxent)
disp mode

stats(fractiles) dtn
scatter(footer="",style=line,hgrid=||%mean,%median,mode||) 1
# xs mu 1 ndraws 1

Re 4. Into the procedure for the LOG and SQRT respectively, I am passing in the transformed actual, and fmean (forecast) & fsigma (forecast std) on the transformed scale.

For the SQRT, I'm not certain of this formula?

Code: Select all

* NON-CENTRAL CHI-SQUARED
set scores startl endl = %logdensity(fsigma**2.0,actual-fcast)-log(2)-0.5*log(actual)

[TomDoan]

I still don't understand why you would overwrite MU with the density estimates. The two have nothing to do with each other. MU has a purpose. Use a different series name for the densities. (FS maybe? Have you looked at the examples of the use of DENSITY?)

You have a Gaussian model for sqrt(actual). That's what you need to use inside the %LOGDENSITY, not actual by itself.

[ac_1]

Sorry I was rushing.

Code: Select all

* empirical distribution

clear memory
*seed 71271

comp ndraws = 1000
set mu 1 ndraws = 5.0
set sigma 1 ndraws = 0.5
set dtn 1 ndraws = (mu+sigma*%ran(1.0))

density(smoothing=2.0) dtn 1 ndraws xs fs

prin / fs

extremum fs

comp mode = xs(%maxent)
disp 'mode' mode

stats(fractiles) dtn
scatter(footer="",style=line,hgrid=||%mean,%median,mode||) 1
# xs fs 1 ndraws 1

Regarding LPDS

(i) NORMAL LPDS I get reasonable LPDS results!

(ii) & (iii) For both

Code: Select all

* LOG-NORMAL
set scores startl endl = ( 1.0/(actual*fsigma*sqrt(2.0*%pi)) * exp(-((log(actual)-fcast)^2.0)/(2.0*fsigma^2.0)) )
set lscores startl endl = log(scores)

* NON-CENTRAL CHI-SQUARED
set scores startl endl = %logdensity(fsigma**2.0,sqrt(actual)-fcast)-log(2)-0.5*log(actual)
set lscores startl endl = log(scores)

sometimes I get negative scores and I cannot calculate lscores (the log of a negative value), therefore I get 0 for LPDS - which is why I am questioning the scores formulae.

I am passing in the transformed actual, and fmean (forecast) & fsigma (forecast std) on the transformed scale.

Attached the procedure for LPDS.

[TomDoan]

No. No. No. The statistic is the log density. Not the log of the log density. Everything would be simpler if you would use a consistent scheme for evaluating the log densities:

Code: Select all

set(first=0.0) x 1 100 = .9*x{1}+%ran(1.0)
set x = x+6.0
*
boxjenk(const,ar=1,define=lineareq) x 1 80
uforecast(equation=lineareq,static) linearf 81 100
sstats 81 100 %logdensity(%sigmasq,x-linearf)>>linlpds
*
set sqrtx = sqrt(x)
boxjenk(const,ar=1,define=sqrteq) sqrtx 1 80
uforecast(equation=sqrteq,static) sqrtf 81 100
sstats 81 100 %logdensity(%sigmasq,sqrt(x)-sqrtf)-log(2*sqrt(x))>>sqrtlpds
*
set logx = log(x)
boxjenk(const,ar=1,define=logeq) logx 1 80
uforecast(equation=logeq,static) logf 81 100
sstats 81 100 %logdensity(%sigmasq,log(x)-logf)-log(x)>>loglpds
?linlpds sqrtlpds loglpds

[ac_1]

Attached

LPDS___examples_multi-step-ahead.rpf
acLPDS.src

LPDS results are the same using acLPDS and as your examples in a main program (just the change the source path).

The procedure has %SIGMASQ the fitted maximum likelihood variance estimate as a (REAL).

I want to calculate LPDS for recursive/rolling one-step head forecasts. Generating the forecasts is easy. How to handle/aggregate the sigmasq's within the procedure for all 3 examples: linear, log, sqrt?

[TomDoan]

Aren't the variances time-varying, just like the forecasts? Put them in a series and pass it through to your procedure.