Calc

Calc is a Fortran-based interactive statistics interpreter with a session-to-Fortran transpiler.

What the interpreter supports

• Scalars and 1D real arrays.
• Array literals, slicing, arithmetic, reductions, and basic control flow.
• Random simulation, fitting, and properties of many probability distributions.
• Time-series helpers including ACF/PACF and AR/MA/ARMA/ARFIMA utilities.
• Robust summaries, common statistical tests, and AIC-based distribution scanning.
• Plotting via gnuplot.

Some interpreter code examples:

! General
run("code.txt")
calc code.txt

! Scalars and arithmetic
n = 10
r = n / 2
r^3
const m = 10^3

! Vectors
y = [1, 2, 3, 4, 5]
v = 10 * arange(10)
t = grid(11, 0.0, 0.1)

! Element-wise operations
z = n * y
w = y + v(1:size(y))
z ^ 2

! Slicing
v([2 4 6 8])
v(3:9:2)
reverse(v)

! Random numbers and time series
x = runif(10)                                ! 10 iid Uniform(0,1) draws
x0 = runif()                                 ! one Uniform(0,1) draw
rn = rnorm(5)                                ! 5 iid standard Normal draws
tn = rnct(5, 8.0, 1.5)                       ! 5 iid noncentral t draws (df=8, ncp=1.5)
mx = rmixnorm(5, [0.7,0.3], [0.0,2.0], [1.0,0.5]) ! 5 draws from a 2-component normal mixture
arsim(1000, [0.5, -0.4])                     ! simulate AR(2) series
acf(x, 10)                                   ! sample ACF for lags 1..10
acf(x, 10, plot=.true.)                      ! sample ACF and plot
pacf(x, 10)                                  ! sample PACF for lags 1..10
pacf(x, 10, plot=.true.)                     ! sample PACF and plot
acfpacf(x, 10)                               ! print aligned ACF/PACF table
acfpacf(x, 10, plot=.true.)                  ! table + joint ACF/PACF plot
acfpacfar(x, 10)                             ! print ACF/PACF/AR-coefficient table
acfpacfar(x, 10, plot=.true.)                ! table + joint ACF/PACF/AR plot
fiacf(0.25, 10)                              ! theoretical ACF of ARFIMA(0,d,0)
arfimaacf([0.4], [0.2], 0.25, 10)            ! theoretical ACF of ARFIMA(1,d,1)
fracdiff(x, 0.3)                             ! fractional differencing (1-L)^0.3 x
arfimafit(x, 1, 1)                           ! fit ARFIMA(1,d,1)
arfimasim(1000, 0.25, phi=[0.4], theta=[0.2]) ! simulate ARFIMA(1,d,1)
polyroots([-6, -1, 1])                       ! roots of x^2 - x - 6 -> packed as [Re1, Im1, Re2, Im2]
armastab([0.6, 0.2], [0.3])                  ! [is_stationary, is_invertible, min_mod_ar, min_mod_ma]
armastab(ma=[0.3])                           ! MA-only invertibility check
adf_stat(x, 4)                               ! ADF t-statistic with 4 lagged differences
adf(x, 4)                                    ! ADF report with p-value and critical values
phillips_perron_stat(x, 8)                   ! PP tau-statistic with HAC bandwidth 8
phillips_perron(x, 8)                        ! PP report with p-value and critical values

! Stats
sum(x)
product(x)
mean(x)
trimmean(x, 0.1)      ! 10% trimmed mean (drop lowest/highest 10%)
winsor_mean(x, 0.1)   ! 10% winsorized mean (cap tails at 10% quantiles)
mad(x)                ! median absolute deviation (robust scale)
iqr(x)                ! interquartile range
iqr_scale(x)          ! robust sd estimate: IQR / 1.349
geomean(x)
harmean(x)
sd(x)
[mean(x) sd(x) minval(x) maxval(x)]
mssk(x)                             ! mean, sd, skew, kurtosis
mssk_norm(0, 1)                     ! theoretical Normal moments [mean, sd, skew, kurtosis]
median(x)
rank(x)
stdz(x)

! Cumulative and differencing
cumsum(y)
cumprod(y)
diff(y)

! Sorting and ordering
sort(x)
p = indexx(y)                     ! permutation that sorts y
y(p)                              ! y reordered by that permutation (same values as sort(y))

! Head/tail
head(v)                            ! first 5 values (default)
head(v, 3)                         ! first 3 values
tail(v)                            ! last 5 values (default)
tail(v, 3)                         ! last 3 values

! Comparisons
x > 0.5
x <= maxval(x)
y == [1 2 3 4 5]
y /= 3
y >= 4

! Two-vector functions
cor(x, y)                         ! Pearson correlation between x and y
cor(x, y, method=spearman)        ! Spearman rank correlation
cor(x, y, method=kendall)         ! Kendall tau-b correlation
cor(x, y, method=["pearson", "spearman", "kendall"]) ! all three pairwise in one call
cov(x, y)                         ! sample covariance between x and y
cor                               ! labeled correlation matrix for all same-length vectors in workspace
cor(x, y, z)                      ! labeled correlation matrix for the listed vectors
cor(x, y, z, method=kendall)      ! labeled matrix with Kendall method
cor(x, y, z, method=["pearson","kendall"]) ! one labeled matrix per method
dot(x, y)                         ! dot product of x and y
min(x, y)                         ! element-wise minimum of x and y
max(x, 0.5)                       ! element-wise maximum of x and scalar 0.5
ttest2(x, y)                      ! Welch two-sample t test -> [t, df, p]
ks2_test(x, y)                    ! two-sample KS test -> [D, p]

! Workspace
?vars                              ! list currently defined variables and their values
read prices.csv                    ! read columns from prices.csv into workspace variables
clear                              ! remove user-defined variables and user-defined functions

! Control flow
if (mean(x) > 0.5) then
  "high"
else
  "low"
end if

! Do loop
do i=1,5
  i, i^2
end do
do i=1,5 i,i^2                     ! one-line do loop

i = 0
do                                 ! potentially infinite loop; include an exit condition
  i = i + 1
  if (i > 10) exit
end do

! For loop over collection
for z in [0.1, 0.2, 0.3]
  z, sqrt(z)
end for
for z in [0.1, 0.2, 0.3] z,sqrt(z) ! one-line for loop

! User-defined functions (arguments are read-only / intent(in)-style)
function center(a)
  center = a - mean(a)
end function
xc = center(x)

function xnorm(a, power=2)         ! default argument in user function
  xnorm = sum(abs(a)^power)^(1/power)
end function
xnorm(x)                            ! uses power=2
xnorm(x, 1)                         ! overrides default
xnorm(x= x, power=1)                ! named arguments

! User-defined subroutines (default argument intent is inout)
subroutine bump(a, b)
  a = a + 1
  b = 2*b
end subroutine
call bump(x, y)

subroutine sum2(a, b, c)
  intent(in) :: a, b
  intent(out) :: c
  c = a + b
end subroutine
call sum2(2, 3, z)

subroutine shift_scale(y, shift=0, scale=1)
  intent(in out) :: y
  intent(in) :: shift, scale        ! defaults allowed for intent(in)
  y = scale*y + shift
end subroutine
call shift_scale(x)                 ! uses shift=0, scale=1
call shift_scale(x, 1.5, 0.5)       ! overrides defaults
call shift_scale(y=x, shift=1.5)    ! named arguments

acf/pacf return lags 1..n and plotting is optional (plot=.false. by default). acfpacf/acfpacfar can also optionally plot. adf/phillips_perron print unit-root test reports; adf_stat/phillips_perron_stat return just the test statistic. head/tail accept an optional second argument for the number of elements to return. One-line do/for loop bodies must be a single statement (you can still use ; within that statement). Defaults in function/subroutine headers must be trailing. In subroutines, defaults are allowed only for arguments declared intent(in). Named arguments are supported in user-defined function and subroutine calls; positional arguments cannot appear after a named argument.

Regression and model fitting

! Simple linear regression
x = runif(200)
y = 1.0 + 2.0*x + 0.2*rnorm(200)
regress(x, y)                     ! with intercept by default
regress(x, y, intcp=0)            ! no-intercept regression
huber_regress(y, x)               ! robust simple regression (Huber loss)
huber_regress(y, x, c=1.345)      ! robust threshold tuning
bisquare_regress(y, x)            ! robust simple regression (Tukey bisquare)
bisquare_regress(y, x, c=4.685)   ! robust threshold tuning
dist_regress(normal, y, x)        ! Gaussian-error regression
dist_regress(t, y, x)             ! Student-t regression, df estimated if omitted
dist_regress(t, y, x, df=[4,6,8,12]) ! choose best df from candidate grid

! Multiple regression
z = x^2
regress(y, x, z)                  ! multiple predictors via regress(...)
poly1reg(y, x, 3)                 ! polynomial regression in one predictor up to degree 3
dist_regress(normal, y, x, z)     ! Gaussian multiple regression
dist_regress(t, y, x, z, df=8)    ! t-errors multiple regression with fixed df

! No-predictor mode
dist_regress(normal, y)           ! intercept-only normal model
dist_regress(t, y)                ! intercept-only t model (df estimated)
dist_regress(normal, y, intcp=0)  ! zero-mean normal model
dist_regress(t, y, intcp=0, df=8) ! zero-mean t model with fixed df

! AR/MA/ARMA fitting helpers
arfit(y, 1, 5)                    ! fit AR orders 1..5 and report fit metrics
mafit(y, 1, 5)                    ! fit MA orders 1..5 and report fit metrics
armafit(y, 1, 1)                  ! fit one ARMA(1,1) model
arsimfit(10^4, [0.6, 0.4])        ! simulate AR and fit default order=size(phi)
masimfit(10^4, [0.6, 0.4])        ! simulate MA and fit default order=size(theta)
armasimfit(10^4, [0.6], [0.3])    ! simulate ARMA and fit default p=size(ar), q=size(ma)
armafitgrid(y, 0, 3, 0, 3)        ! grid search over ARMA(p,q), p=0..3 and q=0..3
armafitaic(y, 0, 5, 0, 5)         ! choose ARMA order by information criterion over p,q ranges
fit_mixnorm(y, 2)                 ! fit a 2-component normal mixture [wgt, mean, sd]
fit_mixnorm_aic(y, 1, 5, nstart=5, verbose=.true., plot=.true.) ! choose mixture size by AIC

arfimafit(x, p, q) prints a fit table including npar (number of estimated parameters), RMSE/AIC/BIC, and parameter estimates. dist_regress prints distribution, sample size,