{UTF8}
:MACRO: TREND2(x,horizon,smoothing time)
TREND2 = smoothed error/horizon
~ x/horizon
~ Second order trend estimate for rate of change of input x, with 0 initial \
trend

smoothed error = SMOOTH(error,horizon)
~ x
~ 
smoothed x = SMOOTH(x,smoothing time)
~ x
~ 
error = x  smoothed x
~ x
~ 
:END OF MACRO:
:MACRO: INIT(x)
INIT = INITIAL(x)
~ x
~ 
:END OF MACRO:
:MACRO: PINK NOISE(mean, std deviation, correlation time, seed)
PINK NOISE = INTEG(updating pink noise,mean+std deviation*RANDOM NORMAL(6,6,0,1,seed\
))
~ mean
~ Contributed by Ed Anderson, MIT/U. Texas  Austin
Description: The pink noise molecule described generates a simple random
series with autocorrelation. This is useful in representing time series,
like rainfall from day to day, in which today's value has some correlation
with what happened yesterday. This particular formulation will also have
properties such as standard deviation and mean that are insensitive both to
the time step and the correlation (smoothing) time. Finally, the output as
a whole and the difference in values between any two days is guaranteed to
be Gaussian (normal) in distribution.
Behavior: Pink noise series will have both a historical and a random
component during each period. The relative "trendtonoise" ratio is
controlled by the length of the correlation time. As the correlation time
approaches zero, the pink noise output will become more independent of its
historical value and more "noisy." On the other hand, as the correlation
time approaches infinity, the pink noise output will approximate a
continuous time random walk or Brownian motion. Displayed above are two
time series with correlation times of 1 and 8 months. While both series
have approximately the same standard deviation, the 1month correlation time
series is less smooth from period to period than the 8month series, which
is characterized by "sustained" swings in a given direction. Note that this
behavior will be independent of the timestep.
The "pink" in pink noise refers to the power spectrum of the output. A time
series in which each period's observation is independent of the past is
characterized by a flat or "white" power spectrum. Smoothing a time series
attenuates the higher or "bluer" frequencies of the power spectrum, leaving
the lower or "redder" frequencies relatively stronger in the output.
Caveats: This assumes the use of Euler integration with a time step of no
more than 1/4 of the correlation time. Very long correlation times should be
avoided also as the multiplication in the scaled white noise will become
progressively less accurate.
Technical Notes: This particular form of pink noise is superior to that of
Britting presented in Richardson and Pugh (1981) because the Gaussian
(Normal) distribution of the output does not depend on the Central Limit
Theorem. (Dynamo did not have a Gaussian random number generator and hence
R&P had to invoke the CLM to get a normal distribution.) Rather, this
molecule's normal output is a result of the observations being a sum of
Gaussian draws. Hence, the series over short intervals should better
approximate normality than the macro in R&P.
MEAN: This is the desired mean for the pink noise.
STD DEVIATION: This is the desired standard deviation for the pink noise.
CORRELATION TIME: This is the smooth time for the noise, or for the more technically \
minded this is the inverse of the filter's cutoff frequency in radians.
Updated by Tom Fiddaman, 2010, to include a random initial value,
correct units, and use TIME STEP$ keyword

updating pink noise = gap/correlation time
~ mean/correlation time
~ 
gap = scaled white noisePINK NOISE
~ mean
~ 
scaled white noise =mean+white noise*std deviation*SQRT((2time step$/correlation time
)/(time step$/correlation time))
~ mean
~ This adjusts the standard deviation of the white noise to compensate for the time \
step and the
correlation time to produce the appropriate pink noise std \
deviation.

white noise = RANDOM NORMAL(6,6,0,1,seed)
~ dmnl
~ This is an independent, identically distributed random quantity drawn every time \
step. The distribution is gaussian with mean = 0 and variance = 1.
Note that RANDOM NORMAL is truncated +/ 6 standard deviations here.
For Vensim 1.62 syntax, remove the arguments to RANDOM NORMAL.

:END OF MACRO:
Smoothing Time=
5
~ Month
~ 
Horizon=
20
~ Month
~ 
Flow Trend=
TREND2(Measured Flow,Horizon,Smoothing Time)
~ drips/(Month*Month)
~ 
Stock Trend=
TREND2(Measured Stock,Horizon,Smoothing Time)
~ drips/Month
~ 
Lagged Measured Flow=
DELAY FIXED(Measured Flow,TIME STEP*"Lag #",Measured Flow)
~ drips/Month
~ 
"Lag #"=
1
~ dmnl [1,5,1]
~ 
Initial Stock=
0
~ drips
~ 
Fractional Outflow Rate=
0
~ fraction/Month [0,?]
~ 
Outflow=
Stock*Fractional Outflow Rate
~ drips/Month
~ 
Flow Abs Error SD=
0
~ drips/Month [0,?]
~ 
Measured Stock=
Stock*EXP(Stock Frac Measurement Error SD*RANDOM NORMAL(6,6,0,1,0))
+Stock Abs Error SD*RANDOM NORMAL(6,6,0,1,0)
~ drips
~ 
Measured Flow=
Flow*EXP(Flow Frac Measurement Error SD*RANDOM NORMAL(6,6,0,1,0))
+Flow Abs Error SD*RANDOM NORMAL(6,6,0,1,0)
~ drips/Month
~ 
Stock Abs Error SD=
0
~ drips [0,?]
~ 
Stock= INTEG (
Flow+NoiseOutflow,
Initial Stock)
~ drips
~ 
Base Flow=
1
~ drips/Month
~ 
Correlation Time=
1
~ Month [1,100]
~ 
Driving Noise SD=
0
~ drips/Month [0,?]
~ Standard deviation of driving noise.

Flow=
Base Flow + STEP(Flow Step, Step Time) + RAMP( Slope , Ramp Time , FINAL TIME )
+ Random Flow SD*PINK NOISE(0,1,Correlation Time,0)
~ drips/Month
~ 
Flow Frac Measurement Error SD=
0
~ fraction [0,?]
~ 
Flow Step=
0
~ drips/Month
~ 
Noise=
Driving Noise SD*RANDOM NORMAL(6,6,0,1,0)
~ drips/Month
~ 
NOISE SEED=
1
~ dmnl [1,10000,1]
~ 
Ramp Time=
20
~ Month
~ 
Random Flow SD=
0
~ drips/Month [0,1]
~ 
Slope=
0
~ drips/Month/Month
~ 
Step Time=
10
~ Month [0,100]
~ 
Stock 1st Difference=
Measured StockSMOOTH(Measured Stock,TIME STEP)
~ drips
~ 
Stock Frac Measurement Error SD=
0
~ fraction [0,?]
~ 
********************************************************
.Control
********************************************************~
Simulation Control Parameters

FINAL TIME = 100
~ Month
~ The final time for the simulation.

INITIAL TIME = 0
~ Month
~ The initial time for the simulation.

SAVEPER =
TIME STEP
~ Month [0,?]
~ The frequency with which output is stored.

TIME STEP = 1
~ Month [0,?]
~ The time step for the simulation.

\\\/// Sketch information  do not modify anything except names
V300 Do not put anything below this section  it will be ignored
*View 1
$192192192,0,Arial120000000025511111196,96,100,0
10,1,Stock,501,274,40,20,3,3,0,0,0,0,0,0
12,2,48,289,274,10,8,0,3,0,0,1,0,0,0
1,3,5,1,4,0,0,22,1,0,0,000,12000,1(423,274)
1,4,5,2,100,0,0,22,1,0,0,000,12000,1(336,274)
11,5,48,380,274,6,8,34,3,0,0,1,0,0,0
10,6,Flow,380,293,18,11,40,3,0,2,1,0,0,0,000,000,1201280
12,7,48,499,421,10,8,0,3,16,2,1,0,0,0,000,000,121280128
1,8,10,1,4,16,0,22,1,0,0,000,12000,1(497,320)
1,9,10,7,100,16,0,22,1,0,0,000,12000,1(497,386)
11,10,48,497,353,8,6,33,3,16,0,4,0,0,0
10,11,Noise,528,353,23,10,40,3,16,2,1,0,0,0,000,000,121280128
10,12,Measured Stock,503,162,59,10,8,3,13,2,0,0,0,0,000,000,1225500
10,13,Measured Flow,334,197,56,10,8,3,16,2,0,0,0,0,000,000,1225500
1,14,5,13,0,16,0,0,1,64,0,000,12000,1(361,243)
1,15,1,12,0,13,0,0,1,64,0,000,12000,1(501,219)
10,16,Base Flow,182,296,39,10,8,3,12,2,1,0,0,0,000,000,1201280
1,17,16,6,0,12,0,0,1,64,0,000,12000,1(284,294)
10,18,Correlation Time,521,500,41,18,8,3,12,2,1,0,0,0,000,000,1201280
1,19,18,6,0,12,0,0,1,64,0,000,12000,1(451,398)
10,20,FINAL TIME,177,335,55,11,8,2,12,3,1,0,0,0,128128128,000,1201280
1,21,20,6,0,12,0,0,1,64,0,000,12000,1(289,311)
10,22,Random Flow SD,420,458,51,18,8,3,12,2,1,0,0,0,000,000,1201280
1,23,22,6,0,12,0,0,1,64,0,000,12000,1(400,378)
10,24,Flow Step,215,410,37,10,8,3,12,2,1,0,0,0,000,000,1201280
1,25,24,6,0,12,0,0,1,64,0,000,12000,1(290,356)
10,26,Ramp Time,328,485,43,10,8,3,12,2,1,0,0,0,000,000,1201280
1,27,26,6,0,12,0,0,1,64,0,000,12000,1(351,396)
10,28,Slope,254,452,23,10,8,3,12,2,1,0,0,0,000,000,1201280
1,29,28,6,0,12,0,0,1,64,0,000,12000,1(311,378)
10,30,Step Time,191,370,38,10,8,3,12,2,1,0,0,0,000,000,1201280
1,31,30,6,0,12,0,0,1,64,0,000,12000,1(282,333)
12,32,0,999,182,150,150,3,44,0,0,1,0,0,0
Stock_vs_Flow
10,33,Driving Noise SD,599,442,50,18,8,3,16,2,1,0,0,0,000,000,121280128
1,34,33,11,0,16,0,0,1,64,0,000,12000,1(564,398)
10,35,NOISE SEED,750,499,51,10,8,3,16,0,0,0,0,0
10,36,Stock Frac Measurement Error SD,723,89,71,27,8,3,16,2,0,0,0,0,000,000,1225500
1,37,36,12,0,16,0,0,1,64,0,000,12000,1(599,129)
10,38,Flow Frac Measurement Error SD,142,188,71,27,8,3,13,2,1,0,0,0,000,000,1225500
1,39,38,13,0,16,0,0,1,64,0,000,12000,1(238,192)
10,40,Stock 1st Difference,746,164,39,18,8,3,16,2,0,0,0,0,000,000,122551280
1,41,12,40,0,16,0,0,1,64,0,000,12000,1(627,162)
10,42,TIME STEP,746,224,50,11,8,2,16,3,1,0,0,0,128128128,000,12128128128
1,43,42,40,0,16,0,0,1,64,0,000,12000,1(746,204)
12,44,0,999,481,150,150,3,44,0,0,1,0,0,0
Stock_1st_Difference
10,45,Lagged Measured Flow,239,126,57,18,8,3,13,2,0,0,0,0,000,000,122551280
1,46,13,45,0,13,0,0,1,64,0,000,12000,1(297,169)
10,47,TIME STEP,103,126,50,11,8,2,13,3,1,0,0,0,128128128,000,12128128128
1,48,47,45,0,13,0,0,1,64,0,000,12000,1(160,126)
10,49,Flow Abs Error SD,193,246,35,18,8,3,16,2,0,0,0,0,000,000,1225500
10,50,Stock Abs Error SD,590,81,39,18,8,3,16,2,0,0,0,0,000,000,1225500
1,51,49,13,0,16,0,0,1,64,0,000,12000,1(259,222)
1,52,50,12,0,16,0,0,1,64,0,000,12000,1(547,120)
12,53,48,698,273,10,8,0,3,8,0,1,0,0,0
1,54,56,53,4,8,0,22,1,0,0,000,12000,1(654,273)
1,55,56,1,100,8,0,22,1,0,0,000,12000,1(574,273)
11,56,48,614,273,6,8,34,3,8,0,1,0,0,0
10,57,Outflow,614,292,28,10,40,3,8,0,1,0,0,0
10,58,Fractional Outflow Rate,711,352,47,18,8,3,8,0,0,0,0,0
1,59,58,57,0,8,0,0,1,64,0,000,12000,1(661,321)
1,60,1,57,1,8,0,0,1,64,0,000,12000,1(563,314)
10,61,Initial Stock,580,222,42,10,8,3,0,0,1,0,0,0
1,62,61,1,0,0,0,0,1,64,1,000,12000,1(553,239)
10,63,"Lag #",104,78,22,10,8,3,13,2,0,0,0,0,111,000,122551280
1,64,63,45,0,13,0,0,1,64,0,000,12000,1(150,94)
10,65,Flow Trend,333,81,41,10,8,3,13,2,0,0,0,0,000,000,1200255
10,66,Stock Trend,462,82,44,10,8,3,13,2,0,0,0,0,000,000,1200255
1,67,13,65,0,13,0,0,1,64,0,000,12000,1(333,145)
1,68,12,66,0,13,0,0,1,64,0,000,12000,1(485,128)
10,69,Horizon,400,120,29,10,8,3,13,2,0,0,0,0,000,000,1200255
1,70,69,65,0,13,0,0,1,64,0,000,12000,1(372,103)
1,71,69,66,0,13,0,0,1,64,0,000,12000,1(424,104)
10,72,Smoothing Time,392,40,59,10,8,3,0,2,0,0,0,0,000,000,000255
1,73,72,65,0,0,0,0,1,64,0,000,0000,1(368,56)
1,74,72,66,0,0,0,0,1,64,0,000,0000,1(420,57)
12,75,0,131,581,79,18,8,7,0,0,1,0,0,0
Tom Fiddaman, 2012  in the public domain
12,76,0,363,579,98,15,8,135,0,18,1,0,253,253,111,000,0U00255
http://models.metasd.comhttp://models.metasd.com
///\\\
:GRAPH Stock_vs_Flow
:TITLE Stock vs Flow
:XAXIS Measured Flow
:DOTS
:SCALE
:VAR Measured Stock
:GRAPH Stock_1st_Difference
:TITLE Stock 1st Difference
:XAXIS Lagged Measured Flow
:DOTS
:SCALE
:VAR Stock 1st Difference
:GRAPH Trend_Comparison
:TITLE Trend Comparison
:XDIV 4
:SCALE
:VAR Flow Trend
:SCALE
:VAR Stock Trend
:L<%^E!@
1:RandomFlowNoiseErr.vdf
9:RandomFlowNoiseErr
15:0,0,0,0,0,0
19:100,0
27:2,
34:0,
4:Time
5:Flow Trend
35:Date
36:YYYYMMDD
37:2000
38:1
39:1
40:2
41:0
24:0
25:100
26:100